00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <ppl-config.h>
00024 #include "MIP_Problem.defs.hh"
00025 #include "globals.defs.hh"
00026 #include "Checked_Number.defs.hh"
00027 #include "Row.defs.hh"
00028 #include "Linear_Expression.defs.hh"
00029 #include "Constraint.defs.hh"
00030 #include "Constraint_System.defs.hh"
00031 #include "Constraint_System.inlines.hh"
00032 #include "Generator.defs.hh"
00033 #include "Scalar_Products.defs.hh"
00034 #include <stdexcept>
00035 #include <deque>
00036 #include <algorithm>
00037 #include <cmath>
00038
00039 #ifndef PPL_NOISY_SIMPLEX
00040 #define PPL_NOISY_SIMPLEX 0
00041 #endif
00042
00043 #ifndef PPL_SIMPLEX_USE_MIP_HEURISTIC
00044 #define PPL_SIMPLEX_USE_MIP_HEURISTIC 1
00045 #endif
00046
00047 #ifdef PPL_NOISY_SIMPLEX
00048 #include <iostream>
00049 #endif
00050
00051
00052 namespace PPL = Parma_Polyhedra_Library;
00053
00054 #if PPL_NOISY_SIMPLEX
00055 namespace {
00056
00057 unsigned long num_iterations = 0;
00058
00059 }
00060 #endif // PPL_NOISY_SIMPLEX
00061
00062 PPL::MIP_Problem::MIP_Problem(const dimension_type dim)
00063 : external_space_dim(dim),
00064 internal_space_dim(0),
00065 tableau(),
00066 working_cost(0, Row::Flags()),
00067 mapping(),
00068 base(),
00069 status(PARTIALLY_SATISFIABLE),
00070 pricing(PRICING_STEEPEST_EDGE_FLOAT),
00071 initialized(false),
00072 input_cs(),
00073 first_pending_constraint(0),
00074 input_obj_function(),
00075 opt_mode(MAXIMIZATION),
00076 last_generator(point()),
00077 i_variables() {
00078
00079 if (dim > max_space_dimension())
00080 throw std::length_error("PPL::MIP_Problem::MIP_Problem(dim, cs, obj, "
00081 "mode):\n"
00082 "dim exceeds the maximum allowed "
00083 "space dimension.");
00084 assert(OK());
00085 }
00086
00087 PPL::MIP_Problem::MIP_Problem(const dimension_type dim,
00088 const Constraint_System& cs,
00089 const Linear_Expression& obj,
00090 const Optimization_Mode mode)
00091 : external_space_dim(dim),
00092 internal_space_dim(0),
00093 tableau(),
00094 working_cost(0, Row::Flags()),
00095 mapping(),
00096 base(),
00097 status(PARTIALLY_SATISFIABLE),
00098 pricing(PRICING_STEEPEST_EDGE_FLOAT),
00099 initialized(false),
00100 input_cs(),
00101 first_pending_constraint(0),
00102 input_obj_function(obj),
00103 opt_mode(mode),
00104 last_generator(point()),
00105 i_variables() {
00106
00107 if (dim > max_space_dimension())
00108 throw std::length_error("PPL::MIP_Problem::MIP_Problem(dim, cs, obj, "
00109 "mode):\n"
00110 "dim exceeds the maximum allowed"
00111 "space dimension.");
00112
00113 if (obj.space_dimension() > dim) {
00114 std::ostringstream s;
00115 s << "PPL::MIP_Problem::MIP_Problem(dim, cs, obj,"
00116 << " mode):\n"
00117 << "obj.space_dimension() == "<< obj.space_dimension()
00118 << " exceeds dim == "<< dim << ".";
00119 throw std::invalid_argument(s.str());
00120 }
00121
00122 if (cs.space_dimension() > dim) {
00123 std::ostringstream s;
00124 s << "PPL::MIP_Problem::MIP_Problem(dim, cs, obj, mode):\n"
00125 << "cs.space_dimension == " << cs.space_dimension() << " exceeds dim == "
00126 << dim << ".";
00127 throw std::invalid_argument(s.str());
00128 }
00129 if (cs.has_strict_inequalities())
00130 throw std::invalid_argument("PPL::MIP_Problem::"
00131 "MIP_Problem(d, cs, obj, m):\n"
00132 "cs contains strict inequalities.");
00133
00134 input_cs.insert(input_cs.end(), cs.begin(), cs.end());
00135 assert(OK());
00136 }
00137
00138 void
00139 PPL::MIP_Problem::add_constraint(const Constraint& c) {
00140 if (space_dimension() < c.space_dimension()) {
00141 std::ostringstream s;
00142 s << "PPL::MIP_Problem::add_constraint(c):\n"
00143 << "c.space_dimension() == "<< c.space_dimension() << " exceeds "
00144 "this->space_dimension == " << space_dimension() << ".";
00145 throw std::invalid_argument(s.str());
00146 }
00147 if (c.is_strict_inequality())
00148 throw std::invalid_argument("PPL::MIP_Problem::add_constraint(c):\n"
00149 "c is a strict inequality.");
00150 input_cs.push_back(c);
00151 if (status != UNSATISFIABLE)
00152 status = PARTIALLY_SATISFIABLE;
00153 assert(OK());
00154 }
00155
00156 void
00157 PPL::MIP_Problem::add_constraints(const Constraint_System& cs) {
00158 if (space_dimension() < cs.space_dimension()) {
00159 std::ostringstream s;
00160 s << "PPL::MIP_Problem::add_constraints(cs):\n"
00161 << "cs.space_dimension() == " << cs.space_dimension()
00162 << " exceeds this->space_dimension() == " << this->space_dimension()
00163 << ".";
00164 throw std::invalid_argument(s.str());
00165 }
00166 if (cs.has_strict_inequalities())
00167 throw std::invalid_argument("PPL::MIP_Problem::add_constraints(cs):\n"
00168 "cs contains strict inequalities.");
00169 input_cs.insert(input_cs.end(), cs.begin(), cs.end());
00170 if (status != UNSATISFIABLE)
00171 status = PARTIALLY_SATISFIABLE;
00172 assert(OK());
00173 }
00174
00175 void
00176 PPL::MIP_Problem::set_objective_function(const Linear_Expression& obj) {
00177 if (space_dimension() < obj.space_dimension()) {
00178 std::ostringstream s;
00179 s << "PPL::MIP_Problem::set_objective_function(obj):\n"
00180 << "obj.space_dimension() == " << obj.space_dimension()
00181 << " exceeds this->space_dimension == " << space_dimension()
00182 << ".";
00183 throw std::invalid_argument(s.str());
00184 }
00185 input_obj_function = obj;
00186 if (status == UNBOUNDED || status == OPTIMIZED)
00187 status = SATISFIABLE;
00188 assert(OK());
00189 }
00190
00191 const PPL::Generator&
00192 PPL::MIP_Problem::feasible_point() const {
00193 if (is_satisfiable())
00194 return last_generator;
00195 else
00196 throw std::domain_error("PPL::MIP_Problem::feasible_point():\n"
00197 "*this is not satisfiable.");
00198 }
00199
00200 const PPL::Generator&
00201 PPL::MIP_Problem::optimizing_point() const {
00202 if (solve() == OPTIMIZED_MIP_PROBLEM)
00203 return last_generator;
00204 else
00205 throw std::domain_error("PPL::MIP_Problem::optimizing_point():\n"
00206 "*this doesn't have an optimizing point.");
00207 }
00208
00209 bool
00210 PPL::MIP_Problem::is_satisfiable() const {
00211
00212 switch (status) {
00213 case UNSATISFIABLE:
00214 assert(OK());
00215 return false;
00216 case SATISFIABLE:
00217
00218 case UNBOUNDED:
00219
00220 case OPTIMIZED:
00221 assert(OK());
00222 return true;
00223 case PARTIALLY_SATISFIABLE:
00224 {
00225 if (i_variables.empty())
00226 return is_lp_satisfiable();
00227
00228 const Variables_Set this_variables_set = integer_space_dimensions();
00229 MIP_Problem& x = const_cast<MIP_Problem&>(*this);
00230 Generator p = point();
00231
00232
00233 x.i_variables.clear();
00234 x.is_lp_satisfiable();
00235 if (is_mip_satisfiable(x, p, this_variables_set)) {
00236 x.last_generator = p;
00237 x.status = SATISFIABLE;
00238
00239 x.i_variables = this_variables_set;
00240 return true;
00241 }
00242 else {
00243 x.status = UNSATISFIABLE;
00244
00245 x.i_variables = this_variables_set;
00246 return false;
00247 }
00248 }
00249 }
00250
00251 throw std::runtime_error("PPL internal error");
00252 }
00253
00254 PPL::MIP_Problem_Status
00255 PPL::MIP_Problem::solve() const{
00256 switch (status) {
00257 case UNSATISFIABLE:
00258 assert(OK());
00259 return UNFEASIBLE_MIP_PROBLEM;
00260 case UNBOUNDED:
00261 assert(OK());
00262 return UNBOUNDED_MIP_PROBLEM;
00263 case OPTIMIZED:
00264 assert(OK());
00265 return OPTIMIZED_MIP_PROBLEM;
00266 case SATISFIABLE:
00267
00268 case PARTIALLY_SATISFIABLE:
00269 {
00270 MIP_Problem& x = const_cast<MIP_Problem&>(*this);
00271 if (i_variables.empty()) {
00272
00273 if (is_lp_satisfiable()) {
00274 x.second_phase();
00275 if (x.status == UNBOUNDED)
00276 return UNBOUNDED_MIP_PROBLEM;
00277 else {
00278 assert(x.status == OPTIMIZED);
00279 return OPTIMIZED_MIP_PROBLEM;
00280 }
00281 }
00282 return UNFEASIBLE_MIP_PROBLEM;
00283 }
00284
00285
00286
00287 const Variables_Set this_variables_set = integer_space_dimensions();
00288 x.i_variables.clear();
00289 if (x.is_lp_satisfiable())
00290 x.second_phase();
00291 else {
00292 x.status = UNSATISFIABLE;
00293
00294 x.i_variables = this_variables_set;
00295 return UNFEASIBLE_MIP_PROBLEM;
00296 }
00297 DIRTY_TEMP0(mpq_class, incumbent_solution);
00298 Generator g = point();
00299 bool have_incumbent_solution = false;
00300
00301 MIP_Problem mip_copy(*this);
00302
00303
00304 mip_copy.i_variables.clear();
00305 MIP_Problem_Status mip_status = solve_mip(have_incumbent_solution,
00306 incumbent_solution, g,
00307 mip_copy,
00308 this_variables_set);
00309
00310 x.i_variables = this_variables_set;
00311 switch (mip_status) {
00312 case UNFEASIBLE_MIP_PROBLEM:
00313 x.status = UNSATISFIABLE;
00314 break;
00315 case UNBOUNDED_MIP_PROBLEM:
00316 x.status = UNBOUNDED;
00317
00318
00319 x.last_generator = g;
00320 break;
00321 case OPTIMIZED_MIP_PROBLEM:
00322 x.status = OPTIMIZED;
00323
00324 x.last_generator = g;
00325 break;
00326 }
00327 assert(OK());
00328 return mip_status;
00329 }
00330 }
00331
00332 throw std::runtime_error("PPL internal error");
00333 }
00334
00335 void
00336 PPL::MIP_Problem::add_space_dimensions_and_embed(const dimension_type m) {
00337
00338
00339 if (m > max_space_dimension() - space_dimension())
00340 throw std::length_error("PPL::MIP_Problem::"
00341 "add_space_dimensions_and_embed(m):\n"
00342 "adding m new space dimensions exceeds "
00343 "the maximum allowed space dimension.");
00344 external_space_dim += m;
00345 if (status != UNSATISFIABLE)
00346 status = PARTIALLY_SATISFIABLE;
00347 assert(OK());
00348 }
00349
00350 void
00351 PPL::MIP_Problem
00352 ::add_to_integer_space_dimensions(const Variables_Set& i_vars) {
00353 if (i_vars.space_dimension() > external_space_dim)
00354 throw std::invalid_argument("PPL::MIP_Problem::"
00355 "add_to_integer_space_dimension(i_vars):\n"
00356 "*this and i_vars are dimension"
00357 "incompatible.");
00358 const dimension_type original_size = i_variables.size();
00359 i_variables.insert(i_vars.begin(), i_vars.end());
00360
00361
00362 if (i_variables.size() != original_size && status != UNSATISFIABLE)
00363 status = PARTIALLY_SATISFIABLE;
00364 }
00365
00366 bool
00367 PPL::MIP_Problem::is_in_base(const dimension_type var_index,
00368 dimension_type& row_index) const {
00369 for (row_index = base.size(); row_index-- > 0; )
00370 if (base[row_index] == var_index)
00371 return true;
00372 return false;
00373 }
00374
00375 void
00376 PPL::MIP_Problem::merge_split_variables(dimension_type var_index,
00377 std::vector<dimension_type>&
00378 unfeasible_tableau_rows) {
00379 const dimension_type tableau_nrows = tableau.num_rows();
00380 const dimension_type column = mapping[var_index].second;
00381
00382 for (dimension_type i = 0; i < tableau_nrows; ++i) {
00383
00384
00385 if (base[i] == mapping[var_index].second) {
00386
00387
00388
00389
00390 #ifndef NDEBUG
00391 for (dimension_type j = 0; j < tableau_nrows; ++j) {
00392 dimension_type dummy = 0;
00393 assert(!is_in_base(mapping[var_index].first, dummy));
00394 }
00395 #endif
00396
00397
00398 base[i] = 0;
00399 unfeasible_tableau_rows.push_back(i);
00400 }
00401 }
00402
00403 const dimension_type tableau_cols = tableau.num_columns();
00404
00405 if (column != tableau_cols - 1) {
00406 std::vector<dimension_type> cycle;
00407 for (dimension_type j = tableau_cols - 1; j >= column; --j)
00408 cycle.push_back(j);
00409 cycle.push_back(0);
00410 tableau.permute_columns(cycle);
00411 }
00412 tableau.remove_trailing_columns(1);
00413
00414
00415 mapping[var_index].second = 0;
00416
00417
00418 const dimension_type base_size = base.size();
00419 for (dimension_type i = base_size; i-- > 0; )
00420 if (base[i] > column)
00421 --base[i];
00422 const dimension_type mapping_size = mapping.size();
00423 for (dimension_type i = mapping_size; i-- > 0; ) {
00424 if (mapping[i].first > column)
00425 --mapping[i].first;
00426 if (mapping[i].second > column)
00427 --mapping[i].second;
00428 }
00429 }
00430
00431 bool
00432 PPL::MIP_Problem::is_satisfied(const Constraint& c, const Generator& g) {
00433
00434
00435 int sp_sign = g.space_dimension() <= c.space_dimension()
00436 ? Scalar_Products::sign(g, c)
00437 : Scalar_Products::sign(c, g);
00438 return c.is_inequality() ? sp_sign >= 0 : sp_sign == 0;
00439 }
00440
00441 bool
00442 PPL::MIP_Problem::is_saturated(const Constraint& c, const Generator& g) {
00443
00444
00445 int sp_sign = g.space_dimension() <= c.space_dimension()
00446 ? Scalar_Products::sign(g, c)
00447 : Scalar_Products::sign(c, g);
00448 return sp_sign == 0;
00449 }
00450
00451 bool
00452 PPL::MIP_Problem::parse_constraints(dimension_type& tableau_num_rows,
00453 dimension_type& num_slack_variables,
00454 std::deque<bool>& is_tableau_constraint,
00455 std::deque<bool>& nonnegative_variable,
00456 std::vector<dimension_type>&
00457 unfeasible_tableau_rows,
00458 std::deque<bool>& satisfied_ineqs) {
00459 satisfied_ineqs.clear();
00460 satisfied_ineqs.insert(satisfied_ineqs.end(), input_cs.size(),
00461 false);
00462
00463 const dimension_type cs_num_rows = input_cs.size();
00464 const dimension_type cs_space_dim = external_space_dim;
00465
00466
00467
00468
00469
00470
00471
00472
00473 tableau_num_rows = cs_num_rows;
00474 dimension_type tableau_num_cols = 2*cs_space_dim;
00475 num_slack_variables = 0;
00476
00477
00478
00479
00480 is_tableau_constraint = std::deque<bool> (cs_num_rows, true);
00481
00482
00483
00484 nonnegative_variable = std::deque<bool> (cs_space_dim, false);
00485
00486
00487
00488 const dimension_type mapping_size = mapping.size();
00489 for (dimension_type i = std::min(mapping_size, cs_space_dim+1); i-- > 1; )
00490 if (mapping[i].second == 0) {
00491 nonnegative_variable[i-1] = true;
00492 --tableau_num_cols;
00493 }
00494
00495
00496 for (dimension_type i = cs_num_rows; i-- > first_pending_constraint; ) {
00497 const Constraint& cs_i = input_cs[i];
00498 bool found_a_nonzero_coeff = false;
00499 bool found_many_nonzero_coeffs = false;
00500 dimension_type nonzero_coeff_column_index = 0;
00501 for (dimension_type sd = cs_i.space_dimension(); sd-- > 0; ) {
00502 if (cs_i.coefficient(Variable(sd)) != 0) {
00503 if (found_a_nonzero_coeff) {
00504 found_many_nonzero_coeffs = true;
00505 if (cs_i.is_inequality())
00506 ++num_slack_variables;
00507 break;
00508 }
00509 else {
00510 nonzero_coeff_column_index = sd + 1;
00511 found_a_nonzero_coeff = true;
00512 }
00513 }
00514 }
00515
00516
00517 if (found_many_nonzero_coeffs) {
00518
00519
00520
00521
00522
00523 if (cs_i.is_inequality() && is_satisfied(cs_i, last_generator))
00524 satisfied_ineqs[i] = true;
00525 continue;
00526 }
00527
00528 if (!found_a_nonzero_coeff) {
00529
00530
00531 if (cs_i.is_inequality()) {
00532 if (cs_i.inhomogeneous_term() < 0)
00533
00534 return false;
00535 }
00536 else
00537
00538 if (cs_i.inhomogeneous_term() != 0)
00539
00540 return false;
00541
00542 is_tableau_constraint[i] = false;
00543 --tableau_num_rows;
00544 continue;
00545 }
00546 else {
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 const dimension_type nonzero_var_index = nonzero_coeff_column_index - 1;
00580
00581 const int sgn_a
00582 = sgn(cs_i.coefficient(Variable(nonzero_coeff_column_index-1)));
00583 const int sgn_b = sgn(cs_i.inhomogeneous_term());
00584
00585 if (sgn_a == sgn_b) {
00586 if (cs_i.is_inequality())
00587 ++num_slack_variables;
00588 }
00589
00590 else if (cs_i.is_equality()) {
00591 if (!nonnegative_variable[nonzero_var_index]) {
00592 nonnegative_variable[nonzero_var_index] = true;
00593 --tableau_num_cols;
00594 }
00595 }
00596
00597 else if (sgn_b < 0) {
00598 if (!nonnegative_variable[nonzero_var_index]) {
00599 nonnegative_variable[nonzero_var_index] = true;
00600 --tableau_num_cols;
00601 }
00602 ++num_slack_variables;
00603 }
00604
00605 else if (sgn_a > 0) {
00606
00607
00608 if (!nonnegative_variable[nonzero_var_index]) {
00609 nonnegative_variable[nonzero_var_index] = true;
00610 --tableau_num_cols;
00611 if (nonzero_coeff_column_index < mapping_size)
00612 merge_split_variables(nonzero_coeff_column_index,
00613 unfeasible_tableau_rows);
00614 is_tableau_constraint[i] = false;
00615 }
00616 else
00617 is_tableau_constraint[i] = false;
00618 --tableau_num_rows;
00619 }
00620
00621 else
00622 ++num_slack_variables;
00623 }
00624 }
00625 return true;
00626 }
00627
00628 bool
00629 PPL::MIP_Problem::process_pending_constraints() {
00630 const dimension_type num_original_rows = tableau.num_rows();
00631 dimension_type new_rows = 0;
00632 dimension_type new_slacks = 0;
00633 dimension_type new_var_columns = 0;
00634 std::deque<bool> is_tableau_constraint;
00635 std::deque<bool> nonnegative_variable;
00636 std::vector<dimension_type> unfeasible_tableau_rows;
00637 std::deque<bool> satisfied_ineqs;
00638
00639
00640
00641 if (!parse_constraints(new_rows, new_slacks, is_tableau_constraint,
00642 nonnegative_variable, unfeasible_tableau_rows,
00643 satisfied_ineqs)) {
00644 status = UNSATISFIABLE;
00645 return false;
00646 };
00647
00648 const dimension_type first_free_tableau_index = tableau.num_columns()-1;
00649
00650 if (external_space_dim > internal_space_dim) {
00651 const dimension_type space_diff = external_space_dim - internal_space_dim;
00652 for (dimension_type i = 0, j = 0; i < space_diff; ++i, ++j) {
00653
00654
00655
00656 if (!nonnegative_variable[internal_space_dim+i]) {
00657 mapping.push_back(std::make_pair(first_free_tableau_index+j,
00658 first_free_tableau_index+1+j));
00659 ++j;
00660 new_var_columns += 2;
00661 }
00662
00663 else {
00664 mapping.push_back(std::make_pair(first_free_tableau_index+j, 0));
00665 ++new_var_columns;
00666 }
00667 }
00668 }
00669
00670
00671
00672 dimension_type num_satisfied_ineqs = std::count(satisfied_ineqs.begin(),
00673 satisfied_ineqs.end(),
00674 true);
00675 const dimension_type unfeasible_tableau_rows_size
00676 = unfeasible_tableau_rows.size();
00677 const dimension_type artificial_cols
00678 = new_rows + unfeasible_tableau_rows_size - num_satisfied_ineqs;
00679 const dimension_type new_total_columns
00680 = new_var_columns + new_slacks + artificial_cols;
00681 if (new_rows > 0)
00682 tableau.add_zero_rows(new_rows, Row::Flags());
00683 if (new_total_columns > 0)
00684 tableau.add_zero_columns(new_total_columns);
00685 dimension_type tableau_num_rows = tableau.num_rows();
00686
00687
00688
00689 std::deque<bool> worked_out_row (tableau_num_rows, false);
00690 dimension_type tableau_num_columns = tableau.num_columns();
00691
00692
00693
00694 base.insert(base.end(), new_rows, 0);
00695 const dimension_type base_size = base.size();
00696
00697
00698 dimension_type slack_index = tableau_num_columns - artificial_cols - 1;
00699 dimension_type artificial_index = slack_index;
00700
00701
00702
00703
00704 const dimension_type begin_artificials = artificial_cols > 0
00705 ? artificial_index : 0;
00706
00707
00708 for (dimension_type k = tableau_num_rows, i = input_cs.size();
00709 i-- > first_pending_constraint; )
00710 if (is_tableau_constraint[i]) {
00711
00712 Row& tableau_k = tableau[--k];
00713 const Constraint& cs_i = input_cs[i];
00714 for (dimension_type sd = cs_i.space_dimension(); sd-- > 0; ) {
00715 tableau_k[mapping[sd+1].first] = cs_i.coefficient(Variable(sd));
00716
00717 if (mapping[sd+1].second != 0)
00718 neg_assign(tableau_k[mapping[sd+1].second],
00719 tableau_k[mapping[sd+1].first]);
00720 }
00721 tableau_k[mapping[0].first] = cs_i.inhomogeneous_term();
00722
00723 if (mapping[0].second != 0)
00724 tableau_k[mapping[0].second] = -cs_i.inhomogeneous_term();
00725
00726
00727 if (cs_i.is_inequality()) {
00728 tableau_k[--slack_index] = -1;
00729
00730
00731
00732 if (satisfied_ineqs[i]) {
00733 base[k] = slack_index;
00734 worked_out_row[k] = true;
00735 }
00736 }
00737 for (dimension_type j = base_size; j-- > 0; )
00738 if (k != j && tableau_k[base[j]] != 0 && base[j] != 0)
00739 linear_combine(tableau_k, tableau[j], base[j]);
00740 }
00741
00742
00743
00744
00745 for (dimension_type i = tableau_num_rows; i-- > 0 ; ) {
00746 Row& tableau_i = tableau[i];
00747 if (tableau_i[0] > 0)
00748 for (dimension_type j = tableau_num_columns; j-- > 0; )
00749 neg_assign(tableau_i[j]);
00750 }
00751
00752
00753 working_cost = Row(tableau_num_columns, Row::Flags());
00754
00755
00756 for (dimension_type i = 0; i < unfeasible_tableau_rows_size; ++i) {
00757 tableau[unfeasible_tableau_rows[i]][artificial_index] = 1;
00758 working_cost[artificial_index] = -1;
00759 base[unfeasible_tableau_rows[i]] = artificial_index;
00760 ++artificial_index;
00761 }
00762
00763
00764
00765
00766
00767
00768
00769 for (dimension_type i = num_original_rows; i < tableau_num_rows; ++i) {
00770 if (worked_out_row[i])
00771 continue;
00772 tableau[i][artificial_index] = 1;
00773 working_cost[artificial_index] = -1;
00774 base[i] = artificial_index;
00775 ++artificial_index;
00776 }
00777
00778 const dimension_type end_artificials = artificial_index - 1;
00779
00780
00781
00782 const dimension_type last_obj_index = working_cost.size() - 1;
00783 working_cost[last_obj_index] = 1;
00784
00785
00786 for (dimension_type i = tableau_num_rows; i-- > 0; )
00787 if (working_cost[base[i]] != 0)
00788 linear_combine(working_cost, tableau[i], base[i]);
00789
00790
00791 if (space_dimension() == 0) {
00792 status = OPTIMIZED;
00793 last_generator = point();
00794 return true;
00795 }
00796
00797
00798
00799
00800
00801 if (tableau_num_rows == 0) {
00802 const dimension_type input_obj_function_size
00803 = input_obj_function.space_dimension();
00804 for (dimension_type i = input_obj_function_size; i-- > 0; )
00805
00806
00807
00808
00809 if ((((input_obj_function.coefficient(Variable(i)) > 0
00810 && opt_mode == MAXIMIZATION)
00811 || (input_obj_function.coefficient(Variable(i)) < 0
00812 && opt_mode == MINIMIZATION)) && mapping[i].second == 0)
00813
00814 || (input_obj_function.coefficient(Variable(i)) != 0
00815 && mapping[i].second != 0)) {
00816
00817 last_generator = point(0 * Variable(space_dimension()-1));
00818 status = UNBOUNDED;
00819 return true;
00820 }
00821
00822
00823
00824
00825 status = OPTIMIZED;
00826
00827 last_generator = point(0*Variable(space_dimension()-1));
00828 assert(OK());
00829 return true;
00830 }
00831
00832
00833 bool first_phase_succesful
00834 = (get_control_parameter(PRICING) == PRICING_STEEPEST_EDGE_FLOAT)
00835 ? compute_simplex_using_steepest_edge_float()
00836 : compute_simplex_using_exact_pricing();
00837
00838 #if PPL_NOISY_SIMPLEX
00839 std::cout << "MIP_Problem::solve: 1st phase ended at iteration "
00840 << num_iterations << "." << std::endl;
00841 #endif
00842
00843 if (!first_phase_succesful || working_cost[0] != 0) {
00844
00845 status = UNSATISFIABLE;
00846 return false;
00847 }
00848
00849
00850 if (begin_artificials != 0)
00851 erase_artificials(begin_artificials, end_artificials);
00852 compute_generator();
00853 status = SATISFIABLE;
00854 assert(OK());
00855 return true;
00856 }
00857
00858 namespace {
00859
00860 inline void
00861 assign(double& d, const mpz_class& c) {
00862 d = c.get_d();
00863 }
00864
00865 template <typename T, typename Policy>
00866 inline void
00867 assign(double& d,
00868 const Parma_Polyhedra_Library::Checked_Number<T, Policy>& c) {
00869 d = raw_value(c);
00870 }
00871
00872 }
00873
00874 PPL::dimension_type
00875 PPL::MIP_Problem::steepest_edge_float_entering_index() const {
00876 DIRTY_TEMP0(mpq_class, real_coeff);
00877 const dimension_type tableau_num_rows = tableau.num_rows();
00878 assert(tableau_num_rows == base.size());
00879 double challenger_num = 0.0;
00880 double challenger_den = 0.0;
00881 double current_value = 0.0;
00882 dimension_type entering_index = 0;
00883 const int cost_sign = sgn(working_cost[working_cost.size() - 1]);
00884 for (dimension_type j = tableau.num_columns() - 1; j-- > 1; ) {
00885 const Coefficient& cost_j = working_cost[j];
00886 if (sgn(cost_j) == cost_sign) {
00887
00888
00889 assign(challenger_num, cost_j);
00890 challenger_num = fabs(challenger_num);
00891
00892
00893 challenger_den = 1.0;
00894 for (dimension_type i = tableau_num_rows; i-- > 0; ) {
00895 const Row& tableau_i = tableau[i];
00896 const Coefficient& tableau_ij = tableau_i[j];
00897 if (tableau_ij != 0) {
00898 assert(tableau_i[base[i]] != 0);
00899 assign_r(real_coeff.get_num(), tableau_ij, ROUND_NOT_NEEDED);
00900 assign_r(real_coeff.get_den(), tableau_i[base[i]], ROUND_NOT_NEEDED);
00901 real_coeff.canonicalize();
00902 double float_tableau_value;
00903 assign(float_tableau_value, real_coeff);
00904 challenger_den += float_tableau_value * float_tableau_value;
00905 }
00906 }
00907 double challenger_value = sqrt(challenger_den);
00908
00909
00910 if (entering_index == 0 || challenger_value > current_value) {
00911 current_value = challenger_value;
00912 entering_index = j;
00913 }
00914 }
00915 }
00916 return entering_index;
00917 }
00918
00919 PPL::dimension_type
00920 PPL::MIP_Problem::steepest_edge_exact_entering_index() const {
00921 const dimension_type tableau_num_rows = tableau.num_rows();
00922 assert(tableau_num_rows == base.size());
00923
00924 TEMP_INTEGER(squared_lcm_basis);
00925
00926 std::vector<Coefficient> norm_factor(tableau_num_rows);
00927 {
00928
00929 TEMP_INTEGER(lcm_basis);
00930 lcm_basis = 1;
00931 for (dimension_type i = tableau_num_rows; i-- > 0; )
00932 lcm_assign(lcm_basis, lcm_basis, tableau[i][base[i]]);
00933
00934 for (dimension_type i = tableau_num_rows; i-- > 0; )
00935 exact_div_assign(norm_factor[i], lcm_basis, tableau[i][base[i]]);
00936
00937
00938 lcm_basis *= lcm_basis;
00939 std::swap(squared_lcm_basis, lcm_basis);
00940 }
00941
00942
00943 TEMP_INTEGER(challenger_num);
00944 TEMP_INTEGER(scalar_value);
00945 TEMP_INTEGER(challenger_den);
00946 TEMP_INTEGER(challenger_value);
00947 TEMP_INTEGER(current_value);
00948
00949 TEMP_INTEGER(current_num);
00950 TEMP_INTEGER(current_den);
00951 dimension_type entering_index = 0;
00952 const int cost_sign = sgn(working_cost[working_cost.size() - 1]);
00953 for (dimension_type j = tableau.num_columns() - 1; j-- > 1; ) {
00954 const Coefficient& cost_j = working_cost[j];
00955 if (sgn(cost_j) == cost_sign) {
00956
00957
00958 challenger_num = cost_j * cost_j;
00959
00960
00961 challenger_den = squared_lcm_basis;
00962 for (dimension_type i = tableau_num_rows; i-- > 0; ) {
00963 const Coefficient& tableau_ij = tableau[i][j];
00964
00965 if (tableau_ij != 0) {
00966 scalar_value = tableau_ij * norm_factor[i];
00967 add_mul_assign(challenger_den, scalar_value, scalar_value);
00968 }
00969 }
00970
00971 if (entering_index == 0) {
00972 std::swap(current_num, challenger_num);
00973 std::swap(current_den, challenger_den);
00974 entering_index = j;
00975 continue;
00976 }
00977 challenger_value = challenger_num * current_den;
00978 current_value = current_num * challenger_den;
00979
00980 if (challenger_value > current_value) {
00981 std::swap(current_num, challenger_num);
00982 std::swap(current_den, challenger_den);
00983 entering_index = j;
00984 }
00985 }
00986 }
00987 return entering_index;
00988 }
00989
00990
00991
00992 PPL::dimension_type
00993 PPL::MIP_Problem::textbook_entering_index() const {
00994
00995
00996
00997
00998
00999
01000 const dimension_type cost_sign_index = working_cost.size() - 1;
01001 const int cost_sign = sgn(working_cost[cost_sign_index]);
01002 assert(cost_sign != 0);
01003 for (dimension_type i = 1; i < cost_sign_index; ++i)
01004 if (sgn(working_cost[i]) == cost_sign)
01005 return i;
01006
01007
01008 return 0;
01009 }
01010
01011 void
01012 PPL::MIP_Problem::linear_combine(Row& x,
01013 const Row& y,
01014 const dimension_type k) {
01015 assert(x.size() == y.size());
01016 assert(y[k] != 0 && x[k] != 0);
01017
01018
01019
01020 TEMP_INTEGER(normalized_x_k);
01021 TEMP_INTEGER(normalized_y_k);
01022 normalize2(x[k], y[k], normalized_x_k, normalized_y_k);
01023 for (dimension_type i = x.size(); i-- > 0; )
01024 if (i != k) {
01025 Coefficient& x_i = x[i];
01026 x_i *= normalized_y_k;
01027 #if 1 // CHECKME: the test seems to speed up the GMP computation.
01028 const Coefficient& y_i = y[i];
01029 if (y_i != 0)
01030 sub_mul_assign(x_i, y_i, normalized_x_k);
01031 #else
01032 sub_mul_assign(x_i, y[i], normalized_x_k);
01033 #endif // 1
01034 }
01035 x[k] = 0;
01036 x.normalize();
01037 }
01038
01039
01040 void
01041 PPL::MIP_Problem::pivot(const dimension_type entering_var_index,
01042 const dimension_type exiting_base_index) {
01043 const Row& tableau_out = tableau[exiting_base_index];
01044
01045 for (dimension_type i = tableau.num_rows(); i-- > 0; ) {
01046 Row& tableau_i = tableau[i];
01047 if (i != exiting_base_index && tableau_i[entering_var_index] != 0)
01048 linear_combine(tableau_i, tableau_out, entering_var_index);
01049 }
01050
01051 if (working_cost[entering_var_index] != 0)
01052 linear_combine(working_cost, tableau_out, entering_var_index);
01053
01054 base[exiting_base_index] = entering_var_index;
01055 }
01056
01057
01058 PPL::dimension_type
01059 PPL::MIP_Problem
01060 ::get_exiting_base_index(const dimension_type entering_var_index) const {
01061
01062
01063
01064
01065
01066
01067
01068 const dimension_type tableau_num_rows = tableau.num_rows();
01069 dimension_type exiting_base_index = tableau_num_rows;
01070 for (dimension_type i = 0; i < tableau_num_rows; ++i) {
01071 const Row& t_i = tableau[i];
01072 const int num_sign = sgn(t_i[entering_var_index]);
01073 if (num_sign != 0 && num_sign == sgn(t_i[base[i]])) {
01074 exiting_base_index = i;
01075 break;
01076 }
01077 }
01078
01079 if (exiting_base_index == tableau_num_rows)
01080 return tableau_num_rows;
01081
01082
01083 TEMP_INTEGER(lcm);
01084 TEMP_INTEGER(current_min);
01085 TEMP_INTEGER(challenger);
01086 for (dimension_type i = exiting_base_index + 1; i < tableau_num_rows; ++i) {
01087 const Row& t_i = tableau[i];
01088 const Coefficient& t_ie = t_i[entering_var_index];
01089 const Coefficient& t_ib = t_i[base[i]];
01090 const int t_ie_sign = sgn(t_ie);
01091 if (t_ie_sign != 0 && t_ie_sign == sgn(t_ib)) {
01092 const Row& t_e = tableau[exiting_base_index];
01093 const Coefficient& t_ee = t_e[entering_var_index];
01094 lcm_assign(lcm, t_ee, t_ie);
01095 exact_div_assign(current_min, lcm, t_ee);
01096 current_min *= t_e[0];
01097 abs_assign(current_min);
01098 exact_div_assign(challenger, lcm, t_ie);
01099 challenger *= t_i[0];
01100 abs_assign(challenger);
01101 current_min -= challenger;
01102 const int sign = sgn(current_min);
01103 if (sign > 0
01104 || (sign == 0 && base[i] < base[exiting_base_index]))
01105 exiting_base_index = i;
01106 }
01107 }
01108 return exiting_base_index;
01109 }
01110
01111
01112 bool
01113 PPL::MIP_Problem::compute_simplex_using_steepest_edge_float() {
01114
01115 const unsigned long allowed_non_increasing_loops = 200;
01116 unsigned long non_increased_times = 0;
01117 bool textbook_pricing = false;
01118
01119 TEMP_INTEGER(cost_sgn_coeff);
01120 TEMP_INTEGER(current_num);
01121 TEMP_INTEGER(current_den);
01122 TEMP_INTEGER(challenger);
01123 TEMP_INTEGER(current);
01124
01125 cost_sgn_coeff = working_cost[working_cost.size()-1];
01126 current_num = working_cost[0];
01127 if (cost_sgn_coeff < 0)
01128 neg_assign(current_num);
01129 abs_assign(current_den, cost_sgn_coeff);
01130 assert(tableau.num_columns() == working_cost.size());
01131 const dimension_type tableau_num_rows = tableau.num_rows();
01132
01133 while (true) {
01134
01135 const dimension_type entering_var_index
01136 = textbook_pricing
01137 ? textbook_entering_index()
01138 : steepest_edge_float_entering_index();
01139
01140
01141 if (entering_var_index == 0)
01142 return true;
01143
01144
01145 const dimension_type exiting_base_index
01146 = get_exiting_base_index(entering_var_index);
01147
01148 if (exiting_base_index == tableau_num_rows)
01149 return false;
01150
01151
01152
01153
01154 maybe_abandon();
01155
01156
01157
01158
01159 pivot(entering_var_index, exiting_base_index);
01160
01161
01162
01163 cost_sgn_coeff = working_cost[working_cost.size()-1];
01164
01165 challenger = working_cost[0];
01166 if (cost_sgn_coeff < 0)
01167 neg_assign(challenger);
01168 challenger *= current_den;
01169 abs_assign(current, cost_sgn_coeff);
01170 current *= current_num;
01171 #if PPL_NOISY_SIMPLEX
01172 ++num_iterations;
01173 if (num_iterations % 200 == 0)
01174 std::cout << "Primal Simplex: iteration "
01175 << num_iterations << "." << std::endl;
01176 #endif
01177
01178 assert(challenger >= current);
01179
01180
01181 if (challenger == current) {
01182 ++non_increased_times;
01183
01184
01185 if (non_increased_times > allowed_non_increasing_loops)
01186 textbook_pricing = true;
01187 }
01188
01189
01190 else {
01191 non_increased_times = 0;
01192 textbook_pricing = false;
01193 }
01194 current_num = working_cost[0];
01195 if (cost_sgn_coeff < 0)
01196 neg_assign(current_num);
01197 abs_assign(current_den, cost_sgn_coeff);
01198 }
01199 }
01200
01201 bool
01202 PPL::MIP_Problem::compute_simplex_using_exact_pricing() {
01203 assert(tableau.num_columns() == working_cost.size());
01204 assert(get_control_parameter(PRICING) == PRICING_STEEPEST_EDGE_EXACT
01205 || get_control_parameter(PRICING) == PRICING_TEXTBOOK);
01206
01207 const dimension_type tableau_num_rows = tableau.num_rows();
01208 const bool textbook_pricing
01209 = (PRICING_TEXTBOOK == get_control_parameter(PRICING));
01210
01211 while (true) {
01212
01213 const dimension_type entering_var_index
01214 = textbook_pricing
01215 ? textbook_entering_index()
01216 : steepest_edge_exact_entering_index();
01217
01218 if (entering_var_index == 0)
01219 return true;
01220
01221
01222 const dimension_type exiting_base_index
01223 = get_exiting_base_index(entering_var_index);
01224
01225 if (exiting_base_index == tableau_num_rows)
01226 return false;
01227
01228
01229
01230
01231 maybe_abandon();
01232
01233
01234
01235
01236 pivot(entering_var_index, exiting_base_index);
01237 #if PPL_NOISY_SIMPLEX
01238 ++num_iterations;
01239 if (num_iterations % 200 == 0)
01240 std::cout << "Primal Simplex: iteration "
01241 << num_iterations << "." << std::endl;
01242 #endif
01243 }
01244 }
01245
01246
01247
01248 void
01249 PPL::MIP_Problem::erase_artificials(const dimension_type begin_artificials,
01250 const dimension_type end_artificials) {
01251 const dimension_type tableau_last_index = tableau.num_columns() - 1;
01252 dimension_type tableau_n_rows = tableau.num_rows();
01253
01254 for (dimension_type i = 0; i < tableau_n_rows; ++i)
01255 if (begin_artificials <= base[i] && base[i] <= end_artificials) {
01256
01257 Row& tableau_i = tableau[i];
01258 bool redundant = true;
01259 for (dimension_type j = end_artificials+1; j-- > 1; )
01260 if (!(begin_artificials <= j && j <= end_artificials)
01261 && tableau_i[j] != 0) {
01262 pivot(j, i);
01263 redundant = false;
01264 break;
01265 }
01266 if (redundant) {
01267
01268
01269 --tableau_n_rows;
01270 if (i < tableau_n_rows) {
01271
01272
01273 tableau_i.swap(tableau[tableau_n_rows]);
01274 base[i] = base[tableau_n_rows];
01275 --i;
01276 }
01277 tableau.erase_to_end(tableau_n_rows);
01278 base.pop_back();
01279 }
01280 }
01281
01282
01283
01284
01285
01286 dimension_type num_artificials = 0;
01287 for (dimension_type i = end_artificials + 1; i-- > 1; )
01288 if (begin_artificials <= i && i <= end_artificials) {
01289 ++num_artificials;
01290 tableau.remove_trailing_columns(1);
01291 }
01292
01293
01294 for (dimension_type i = tableau_n_rows; i-- > 0; )
01295 tableau[i][tableau.num_columns()-1] = 0;
01296
01297
01298
01299 working_cost[tableau.num_columns()-1] = working_cost[tableau_last_index];
01300
01301 const dimension_type working_cost_new_size = working_cost.size() -
01302 num_artificials;
01303 working_cost.shrink(working_cost_new_size);
01304 }
01305
01306
01307 void
01308 PPL::MIP_Problem::compute_generator() const {
01309
01310
01311 std::vector<Coefficient> num(external_space_dim);
01312 std::vector<Coefficient> den(external_space_dim);
01313 dimension_type row = 0;
01314
01315 TEMP_INTEGER(lcm);
01316
01317 TEMP_INTEGER(split_num);
01318 TEMP_INTEGER(split_den);
01319
01320
01321 for (dimension_type i = external_space_dim; i-- > 0; ) {
01322 Coefficient& num_i = num[i];
01323 Coefficient& den_i = den[i];
01324
01325
01326 const dimension_type original_var = mapping[i+1].first;
01327 if (is_in_base(original_var, row)) {
01328 const Row& t_row = tableau[row];
01329 if (t_row[original_var] > 0) {
01330 neg_assign(num_i, t_row[0]);
01331 den_i = t_row[original_var];
01332 }
01333 else {
01334 num_i = t_row[0];
01335 neg_assign(den_i, t_row[original_var]);
01336 }
01337 }
01338 else {
01339 num_i = 0;
01340 den_i = 1;
01341 }
01342
01343 const dimension_type split_var = mapping[i+1].second;
01344 if (split_var != 0) {
01345
01346
01347
01348 if (is_in_base(split_var, row)) {
01349 const Row& t_row = tableau[row];
01350 if (t_row[split_var] > 0) {
01351 split_num = -t_row[0];
01352 split_den = t_row[split_var];
01353 }
01354 else {
01355 split_num = t_row[0];
01356 split_den = -t_row[split_var];
01357 }
01358
01359
01360 lcm_assign(lcm, den_i, split_den);
01361 exact_div_assign(den_i, lcm, den_i);
01362 exact_div_assign(split_den, lcm, split_den);
01363 num_i *= den_i;
01364 sub_mul_assign(num_i, split_num, split_den);
01365 if (num_i == 0)
01366 den_i = 1;
01367 else
01368 den_i = lcm;
01369 }
01370
01371
01372 }
01373 }
01374
01375
01376 lcm = den[0];
01377 for (dimension_type i = 1; i < external_space_dim; ++i)
01378 lcm_assign(lcm, lcm, den[i]);
01379
01380
01381 for (dimension_type i = external_space_dim; i-- > 0; ) {
01382 exact_div_assign(den[i], lcm, den[i]);
01383 num[i] *= den[i];
01384 }
01385
01386
01387 Linear_Expression expr;
01388 for (dimension_type i = external_space_dim; i-- > 0; )
01389 expr += num[i] * Variable(i);
01390
01391 MIP_Problem& x = const_cast<MIP_Problem&>(*this);
01392 x.last_generator = point(expr, lcm);
01393 }
01394
01395 void
01396 PPL::MIP_Problem::second_phase() {
01397
01398 assert(status == SATISFIABLE || status == UNBOUNDED || status == OPTIMIZED);
01399
01400 if (status == UNBOUNDED || status == OPTIMIZED)
01401 return;
01402
01403
01404 const dimension_type input_obj_function_sd
01405 = input_obj_function.space_dimension();
01406 Row new_cost(input_obj_function_sd + 1, Row::Flags());
01407 for (dimension_type i = input_obj_function_sd; i-- > 0; )
01408 new_cost[i+1] = input_obj_function.coefficient(Variable(i));
01409 new_cost[0] = input_obj_function.inhomogeneous_term();
01410
01411
01412 if (opt_mode == MINIMIZATION)
01413 for (dimension_type i = new_cost.size(); i-- > 0; )
01414 neg_assign(new_cost[i]);
01415
01416
01417 const dimension_type cost_zero_size = working_cost.size();
01418 Row tmp_cost = Row(new_cost, cost_zero_size, cost_zero_size);
01419 tmp_cost.swap(working_cost);
01420 working_cost[cost_zero_size-1] = 1;
01421
01422
01423 for (dimension_type i = new_cost.size(); i-- > 1; ) {
01424 const dimension_type original_var = mapping[i].first;
01425 const dimension_type split_var = mapping[i].second;
01426 working_cost[original_var] = new_cost[i];
01427 if (mapping[i].second != 0)
01428 working_cost[split_var] = - new_cost[i];
01429 }
01430
01431
01432 for (dimension_type i = tableau.num_rows(); i-- > 0; ) {
01433 const dimension_type base_i = base[i];
01434 if (working_cost[base_i] != 0)
01435 linear_combine(working_cost, tableau[i], base_i);
01436 }
01437
01438 bool second_phase_successful
01439 = (get_control_parameter(PRICING) == PRICING_STEEPEST_EDGE_FLOAT)
01440 ? compute_simplex_using_steepest_edge_float()
01441 : compute_simplex_using_exact_pricing();
01442 compute_generator();
01443 #if PPL_NOISY_SIMPLEX
01444 std::cout << "MIP_Problem::solve: 2nd phase ended at iteration "
01445 << num_iterations << "." << std::endl;
01446 #endif
01447 status = second_phase_successful ? OPTIMIZED : UNBOUNDED;
01448 assert(OK());
01449 }
01450
01451 void
01452 PPL::MIP_Problem
01453 ::evaluate_objective_function(const Generator& evaluating_point,
01454 Coefficient& ext_n,
01455 Coefficient& ext_d) const {
01456 const dimension_type ep_space_dim = evaluating_point.space_dimension();
01457 if (space_dimension() < ep_space_dim)
01458 throw std::invalid_argument("PPL::MIP_Problem::"
01459 "evaluate_objective_function(p, n, d):\n"
01460 "*this and p are dimension incompatible.");
01461 if (!evaluating_point.is_point())
01462 throw std::invalid_argument("PPL::MIP_Problem::"
01463 "evaluate_objective_function(p, n, d):\n"
01464 "p is not a point.");
01465
01466
01467
01468 const dimension_type working_space_dim
01469 = std::min(ep_space_dim, input_obj_function.space_dimension());
01470
01471 const Coefficient& divisor = evaluating_point.divisor();
01472 ext_n = input_obj_function.inhomogeneous_term() * divisor;
01473 for (dimension_type i = working_space_dim; i-- > 0; )
01474 ext_n += evaluating_point.coefficient(Variable(i))
01475 * input_obj_function.coefficient(Variable(i));
01476
01477 normalize2(ext_n, divisor, ext_n, ext_d);
01478 }
01479
01480 bool
01481 PPL::MIP_Problem::is_lp_satisfiable() const {
01482 #if PPL_NOISY_SIMPLEX
01483 num_iterations = 0;
01484 #endif
01485 switch (status) {
01486 case UNSATISFIABLE:
01487 return false;
01488 case SATISFIABLE:
01489
01490 case UNBOUNDED:
01491
01492 case OPTIMIZED:
01493
01494 return true;
01495 case PARTIALLY_SATISFIABLE:
01496 {
01497 MIP_Problem& x = const_cast<MIP_Problem&>(*this);
01498
01499
01500 if (tableau.num_columns() == 0) {
01501
01502
01503 x.tableau.add_zero_columns(2);
01504
01505 x.mapping.push_back(std::make_pair(0, 0));
01506
01507
01508 x.initialized = true;
01509 }
01510
01511
01512 x.process_pending_constraints();
01513
01514 x.first_pending_constraint = input_cs.size();
01515
01516 x.internal_space_dim = x.external_space_dim;
01517 assert(OK());
01518 return (status != UNSATISFIABLE);
01519 }
01520 }
01521
01522 throw std::runtime_error("PPL internal error");
01523 }
01524
01525 PPL::MIP_Problem_Status
01526 PPL::MIP_Problem::solve_mip(bool& have_incumbent_solution,
01527 mpq_class& incumbent_solution_value,
01528 Generator& incumbent_solution_point,
01529 MIP_Problem& lp,
01530 const Variables_Set& i_vars) {
01531
01532 PPL::MIP_Problem_Status lp_status;
01533 if (lp.is_lp_satisfiable()) {
01534 lp.second_phase();
01535 lp_status = (lp.status == OPTIMIZED) ? OPTIMIZED_MIP_PROBLEM
01536 : UNBOUNDED_MIP_PROBLEM;
01537 }
01538 else
01539 return UNFEASIBLE_MIP_PROBLEM;
01540
01541 DIRTY_TEMP0(mpq_class, tmp_rational);
01542
01543 Generator p = point();
01544 TEMP_INTEGER(tmp_coeff1);
01545 TEMP_INTEGER(tmp_coeff2);
01546
01547 if (lp_status == UNBOUNDED_MIP_PROBLEM)
01548 p = lp.last_generator;
01549 else {
01550 assert(lp_status == OPTIMIZED_MIP_PROBLEM);
01551
01552 p = lp.last_generator;
01553 lp.evaluate_objective_function(p, tmp_coeff1, tmp_coeff2);
01554 assign_r(tmp_rational.get_num(), tmp_coeff1, ROUND_NOT_NEEDED);
01555 assign_r(tmp_rational.get_den(), tmp_coeff2, ROUND_NOT_NEEDED);
01556 assert(is_canonical(tmp_rational));
01557 if (have_incumbent_solution
01558 && ((lp.optimization_mode() == MAXIMIZATION
01559 && tmp_rational <= incumbent_solution_value)
01560 || (lp.optimization_mode() == MINIMIZATION
01561 && tmp_rational >= incumbent_solution_value)))
01562
01563 return lp_status;
01564 }
01565
01566 bool found_satisfiable_generator = true;
01567 TEMP_INTEGER(gcd);
01568 const Coefficient& p_divisor = p.divisor();
01569 dimension_type nonint_dim;
01570 for (Variables_Set::const_iterator v_begin = i_vars.begin(),
01571 v_end = i_vars.end(); v_begin != v_end; ++v_begin) {
01572 gcd_assign(gcd, p.coefficient(Variable(*v_begin)), p_divisor);
01573 if (gcd != p_divisor) {
01574 nonint_dim = *v_begin;
01575 found_satisfiable_generator = false;
01576 break;
01577 }
01578 }
01579 if (found_satisfiable_generator) {
01580
01581 if (lp_status == UNBOUNDED_MIP_PROBLEM) {
01582
01583
01584
01585 incumbent_solution_point = p;
01586 return lp_status;
01587 }
01588 if (!have_incumbent_solution
01589 || (lp.optimization_mode() == MAXIMIZATION
01590 && tmp_rational > incumbent_solution_value)
01591 || tmp_rational < incumbent_solution_value) {
01592 incumbent_solution_value = tmp_rational;
01593 incumbent_solution_point = p;
01594 have_incumbent_solution = true;
01595 #if PPL_NOISY_SIMPLEX
01596 TEMP_INTEGER(num);
01597 TEMP_INTEGER(den);
01598 lp.evaluate_objective_function(p, num, den);
01599 std::cerr << "new value found: " << num << "/" << den << std::endl;
01600 #endif
01601 }
01602 return lp_status;
01603 }
01604
01605 assert(nonint_dim < lp.space_dimension());
01606
01607 assign_r(tmp_rational.get_num(), p.coefficient(Variable(nonint_dim)),
01608 ROUND_NOT_NEEDED);
01609 assign_r(tmp_rational.get_den(), p_divisor, ROUND_NOT_NEEDED);
01610 tmp_rational.canonicalize();
01611 assign_r(tmp_coeff1, tmp_rational, ROUND_DOWN);
01612 assign_r(tmp_coeff2, tmp_rational, ROUND_UP);
01613 {
01614 MIP_Problem lp_aux = lp;
01615 lp_aux.add_constraint(Variable(nonint_dim) <= tmp_coeff1);
01616 solve_mip(have_incumbent_solution, incumbent_solution_value,
01617 incumbent_solution_point, lp_aux, i_vars);
01618 }
01619
01620 lp.add_constraint(Variable(nonint_dim) >= tmp_coeff2);
01621 solve_mip(have_incumbent_solution, incumbent_solution_value,
01622 incumbent_solution_point, lp, i_vars);
01623 return have_incumbent_solution ? lp_status : UNFEASIBLE_MIP_PROBLEM;
01624 }
01625
01626 bool
01627 PPL::MIP_Problem::choose_branching_variable(const MIP_Problem& mip,
01628 const Variables_Set& i_vars,
01629 dimension_type& branching_index) {
01630
01631 const Constraint_Sequence& input_cs = mip.input_cs;
01632 const Generator& last_generator = mip.last_generator;
01633 const Coefficient& last_generator_divisor = last_generator.divisor();
01634 Variables_Set candidate_variables;
01635
01636 TEMP_INTEGER(gcd);
01637 for (Variables_Set::const_iterator v_it = i_vars.begin(),
01638 v_end = i_vars.end(); v_it != v_end; ++v_it) {
01639 gcd_assign(gcd,
01640 last_generator.coefficient(Variable(*v_it)),
01641 last_generator_divisor);
01642 if (gcd != last_generator_divisor)
01643 candidate_variables.insert(*v_it);
01644 }
01645
01646 if (candidate_variables.empty())
01647 return true;
01648
01649
01650 const dimension_type input_cs_num_rows = input_cs.size();
01651 std::deque<bool> satisfiable_constraints (input_cs_num_rows, false);
01652 for (dimension_type i = input_cs_num_rows; i-- > 0; )
01653
01654
01655 if (input_cs[i].is_equality()
01656 || is_saturated(input_cs[i], last_generator))
01657 satisfiable_constraints[i] = true;
01658
01659 dimension_type current_num_appearances = 0;
01660 dimension_type winning_num_appearances = 0;
01661
01662
01663
01664 for (Variables_Set::const_iterator v_it = candidate_variables.begin(),
01665 v_end = candidate_variables.end(); v_it != v_end; ++v_it) {
01666 current_num_appearances = 0;
01667 for (dimension_type i = input_cs_num_rows; i-- > 0; )
01668 if (satisfiable_constraints[i]
01669 && *v_it < input_cs[i].space_dimension()
01670 && input_cs[i].coefficient(Variable(*v_it)) != 0)
01671 ++current_num_appearances;
01672 if (current_num_appearances >= winning_num_appearances) {
01673 winning_num_appearances = current_num_appearances;
01674 branching_index = *v_it;
01675 }
01676 }
01677 return false;
01678 }
01679
01680 bool
01681 PPL::MIP_Problem::is_mip_satisfiable(MIP_Problem& lp, Generator& p,
01682 const Variables_Set& i_vars) {
01683
01684 if (!lp.is_lp_satisfiable())
01685 return false;
01686 DIRTY_TEMP0(mpq_class, tmp_rational);
01687
01688 TEMP_INTEGER(tmp_coeff1);
01689 TEMP_INTEGER(tmp_coeff2);
01690 p = lp.last_generator;
01691
01692 bool found_satisfiable_generator = true;
01693 dimension_type nonint_dim;
01694 const Coefficient& p_divisor = p.divisor();
01695
01696 #if PPL_SIMPLEX_USE_MIP_HEURISTIC
01697 found_satisfiable_generator
01698 = choose_branching_variable(lp, i_vars, nonint_dim);
01699 #else
01700 TEMP_INTEGER(gcd);
01701 for (Variables_Set::const_iterator v_begin = i_vars.begin(),
01702 v_end = i_vars.end(); v_begin != v_end; ++v_begin) {
01703 gcd_assign(gcd, p.coefficient(Variable(*v_begin)), p_divisor);
01704 if (gcd != p_divisor) {
01705 nonint_dim = *v_begin;
01706 found_satisfiable_generator = false;
01707 break;
01708 }
01709 }
01710 #endif
01711
01712 if (found_satisfiable_generator)
01713 return true;
01714
01715
01716 assert(nonint_dim < lp.space_dimension());
01717
01718 assign_r(tmp_rational.get_num(), p.coefficient(Variable(nonint_dim)),
01719 ROUND_NOT_NEEDED);
01720 assign_r(tmp_rational.get_den(), p_divisor, ROUND_NOT_NEEDED);
01721 tmp_rational.canonicalize();
01722 assign_r(tmp_coeff1, tmp_rational, ROUND_DOWN);
01723 assign_r(tmp_coeff2, tmp_rational, ROUND_UP);
01724 {
01725 MIP_Problem lp_aux = lp;
01726 lp_aux.add_constraint(Variable(nonint_dim) <= tmp_coeff1);
01727 if (is_mip_satisfiable(lp_aux, p, i_vars))
01728 return true;
01729 }
01730 lp.add_constraint(Variable(nonint_dim) >= tmp_coeff2);
01731 return is_mip_satisfiable(lp, p, i_vars);
01732 }
01733
01734 bool
01735 PPL::MIP_Problem::OK() const {
01736 #ifndef NDEBUG
01737 using std::endl;
01738 using std::cerr;
01739 #endif
01740 const dimension_type input_cs_num_rows = input_cs.size();
01741
01742
01743 for (dimension_type i = input_cs_num_rows; i-- > 0; )
01744 if (!input_cs[i].OK())
01745 return false;
01746
01747 if (!tableau.OK() || !input_obj_function.OK() || !last_generator.OK())
01748 return false;
01749
01750
01751 for (dimension_type i = input_cs_num_rows; i-- > 0; )
01752 if (input_cs[i].is_strict_inequality()) {
01753 #ifndef NDEBUG
01754 cerr << "The feasible region of the MIP_Problem is defined by "
01755 << "a constraint system containing strict inequalities."
01756 << endl;
01757 ascii_dump(cerr);
01758 #endif
01759 return false;
01760 }
01761
01762
01763 if (external_space_dim < input_obj_function.space_dimension()) {
01764 #ifndef NDEBUG
01765 cerr << "The MIP_Problem and the objective function have "
01766 << "incompatible space dimensions ("
01767 << external_space_dim << " < "
01768 << input_obj_function.space_dimension() << ")."
01769 << endl;
01770 ascii_dump(cerr);
01771 #endif
01772 return false;
01773 }
01774
01775 if (status != UNSATISFIABLE && initialized) {
01776
01777
01778 if (external_space_dim != last_generator.space_dimension()) {
01779 #ifndef NDEBUG
01780 cerr << "The MIP_Problem and the cached feasible point have "
01781 << "incompatible space dimensions ("
01782 << external_space_dim << " != "
01783 << last_generator.space_dimension() << ")."
01784 << endl;
01785 ascii_dump(cerr);
01786 #endif
01787 return false;
01788 }
01789
01790 for (dimension_type i = 0; i < first_pending_constraint; ++i)
01791 if (!is_satisfied(input_cs[i], last_generator)) {
01792 #ifndef NDEBUG
01793 cerr << "The cached feasible point does not belong to "
01794 << "the feasible region of the MIP_Problem."
01795 << endl;
01796 ascii_dump(cerr);
01797 #endif
01798 return false;
01799 }
01800
01801
01802
01803 if (!i_variables.empty()) {
01804 TEMP_INTEGER(gcd);
01805 for (Variables_Set::const_iterator v_it = i_variables.begin(),
01806 v_end = i_variables.end(); v_it != v_end; ++v_it) {
01807 gcd_assign(gcd, last_generator.coefficient(Variable(*v_it)),
01808 last_generator.divisor());
01809 if (gcd != last_generator.divisor())
01810 return false;
01811 }
01812 }
01813
01814 const dimension_type tableau_nrows = tableau.num_rows();
01815 const dimension_type tableau_ncols = tableau.num_columns();
01816
01817
01818 if (tableau_nrows != base.size()) {
01819 #ifndef NDEBUG
01820 cerr << "tableau and base have incompatible sizes" << endl;
01821 ascii_dump(cerr);
01822 #endif
01823 return false;
01824 }
01825
01826
01827 if (mapping.size() != external_space_dim + 1) {
01828 #ifndef NDEBUG
01829 cerr << "`input_cs' and `mapping' have incompatible sizes" << endl;
01830 ascii_dump(cerr);
01831 #endif
01832 return false;
01833 }
01834
01835
01836 if (tableau_ncols != working_cost.size()) {
01837 #ifndef NDEBUG
01838 cerr << "tableau and working_cost have incompatible sizes" << endl;
01839 ascii_dump(cerr);
01840 #endif
01841 return false;
01842 }
01843
01844
01845 for (dimension_type i = base.size(); i-- > 0; ) {
01846 if (base[i] > tableau_ncols) {
01847 #ifndef NDEBUG
01848 cerr << "base contains an invalid column index" << endl;
01849 ascii_dump(cerr);
01850 #endif
01851 return false;
01852 }
01853
01854
01855 for (dimension_type j = tableau_nrows; j-- > 0; )
01856 if (i != j && tableau[j][base[i]] != 0) {
01857 #ifndef NDEBUG
01858 cerr << "tableau[i][base[i] must be different from zero" << endl;
01859 ascii_dump(cerr);
01860 #endif
01861 return false;
01862 }
01863 if (tableau[i][base[i]] == 0) {
01864 #ifndef NDEBUG
01865 cerr << "tableau[i][base[j], with i different from j, must not be "
01866 << "a zero" << endl;
01867 ascii_dump(cerr);
01868 #endif
01869 return false;
01870 }
01871 }
01872
01873
01874 for (dimension_type i = tableau_nrows; i-- > 0; )
01875 if (tableau[i][tableau_ncols-1] != 0) {
01876 #ifndef NDEBUG
01877 cerr << "the last column of the tableau must contain only"
01878 "zeroes"<< endl;
01879 ascii_dump(cerr);
01880 #endif
01881 return false;
01882 }
01883 }
01884
01885
01886 return true;
01887 }
01888
01889 void
01890 PPL::MIP_Problem::ascii_dump(std::ostream& s) const {
01891 using namespace IO_Operators;
01892 s << "\nexternal_space_dim: " << external_space_dim << " \n";
01893 s << "\ninternal_space_dim: " << internal_space_dim << " \n";
01894
01895 const dimension_type input_cs_size = input_cs.size();
01896
01897 s << "\ninput_cs( " << input_cs_size << " )\n";
01898 for (dimension_type i = 0; i < input_cs_size; ++i)
01899 input_cs[i].ascii_dump(s);
01900
01901 s << "\nfirst_pending_constraint: " << first_pending_constraint
01902 << std::endl;
01903
01904 s << "\ninput_obj_function\n";
01905 input_obj_function.ascii_dump(s);
01906 s << "\nopt_mode "
01907 << (opt_mode == MAXIMIZATION ? "MAXIMIZATION" : "MINIMIZATION") << "\n";
01908
01909 s << "\ninitialized: " << (initialized ? "YES" : "NO") << "\n";
01910 s << "\npricing: ";
01911 switch (pricing) {
01912 case PRICING_STEEPEST_EDGE_FLOAT:
01913 s << "PRICING_STEEPEST_EDGE_FLOAT";
01914 break;
01915 case PRICING_STEEPEST_EDGE_EXACT:
01916 s << "PRICING_STEEPEST_EDGE_EXACT";
01917 break;
01918 case PRICING_TEXTBOOK:
01919 s << "PRICING_TEXTBOOK";
01920 break;
01921 }
01922 s << "\n";
01923
01924 s << "\nstatus: ";
01925 switch (status) {
01926 case UNSATISFIABLE:
01927 s << "UNSATISFIABLE";
01928 break;
01929 case SATISFIABLE:
01930 s << "SATISFIABLE";
01931 break;
01932 case UNBOUNDED:
01933 s << "UNBOUNDED";
01934 break;
01935 case OPTIMIZED:
01936 s << "OPTIMIZED";
01937 break;
01938 case PARTIALLY_SATISFIABLE:
01939 s << "PARTIALLY_SATISFIABLE";
01940 break;
01941 }
01942 s << "\n";
01943
01944 s << "\ntableau\n";
01945 tableau.ascii_dump(s);
01946 s << "\nworking_cost( " << working_cost.size()<< " )\n";
01947 working_cost.ascii_dump(s);
01948
01949 const dimension_type base_size = base.size();
01950 s << "\nbase( " << base_size << " )\n";
01951 for (dimension_type i = 0; i != base_size; ++i)
01952 s << base[i] << ' ';
01953
01954 s << "\nlast_generator\n";
01955 last_generator.ascii_dump(s);
01956
01957 const dimension_type mapping_size = mapping.size();
01958 s << "\nmapping( " << mapping_size << " )\n";
01959 for (dimension_type i = 1; i < mapping_size; ++i)
01960 s << "\n"<< i << " -> " << mapping[i].first << " -> " << mapping[i].second
01961 << ' ';
01962
01963 s << "\n\ninteger_variables";
01964 i_variables.ascii_dump(s);
01965 }
01966
01967 PPL_OUTPUT_DEFINITIONS(MIP_Problem)
01968
01969 bool
01970 PPL::MIP_Problem::ascii_load(std::istream& s) {
01971 std::string str;
01972 if (!(s >> str) || str != "external_space_dim:")
01973 return false;
01974
01975 if (!(s >> external_space_dim))
01976 return false;
01977
01978 if (!(s >> str) || str != "internal_space_dim:")
01979 return false;
01980
01981 if (!(s >> internal_space_dim))
01982 return false;
01983
01984 if (!(s >> str) || str != "input_cs(")
01985 return false;
01986
01987 dimension_type input_cs_size;
01988
01989 if (!(s >> input_cs_size))
01990 return false;
01991
01992 if (!(s >> str) || str != ")")
01993 return false;
01994
01995 Constraint c(Constraint::zero_dim_positivity());
01996 for (dimension_type i = 0; i < input_cs_size; ++i) {
01997 if (!c.ascii_load(s))
01998 return false;
01999 input_cs.push_back(c);
02000 }
02001
02002 if (!(s >> str) || str != "first_pending_constraint:")
02003 return false;
02004
02005 if (!(s >> first_pending_constraint))
02006 return false;
02007
02008 if (!(s >> str) || str != "input_obj_function")
02009 return false;
02010
02011 if (!input_obj_function.ascii_load(s))
02012 return false;
02013
02014 if (!(s >> str) || str != "opt_mode")
02015 return false;
02016
02017 if (!(s >> str))
02018 return false;
02019
02020 if (str == "MAXIMIZATION")
02021 set_optimization_mode(MAXIMIZATION);
02022 else {
02023 if (str != "MINIMIZATION")
02024 return false;
02025 set_optimization_mode(MINIMIZATION);
02026 }
02027
02028 if (!(s >> str) || str != "initialized:")
02029 return false;
02030 if (!(s >> str))
02031 return false;
02032 if (str == "YES")
02033 initialized = true;
02034 else if (str == "NO")
02035 initialized = false;
02036 else
02037 return false;
02038
02039 if (!(s >> str) || str != "pricing:")
02040 return false;
02041 if (!(s >> str))
02042 return false;
02043 if (str == "PRICING_STEEPEST_EDGE_FLOAT")
02044 pricing = PRICING_STEEPEST_EDGE_FLOAT;
02045 else if (str == "PRICING_STEEPEST_EDGE_EXACT")
02046 pricing = PRICING_STEEPEST_EDGE_EXACT;
02047 else if (str == "PRICING_TEXTBOOK")
02048 pricing = PRICING_TEXTBOOK;
02049 else
02050 return false;
02051
02052 if (!(s >> str) || str != "status:")
02053 return false;
02054
02055 if (!(s >> str))
02056 return false;
02057
02058 if (str == "UNSATISFIABLE")
02059 status = UNSATISFIABLE;
02060 else if (str == "SATISFIABLE")
02061 status = SATISFIABLE;
02062 else if (str == "UNBOUNDED")
02063 status = UNBOUNDED;
02064 else if (str == "OPTIMIZED")
02065 status = OPTIMIZED;
02066 else if (str == "PARTIALLY_SATISFIABLE")
02067 status = PARTIALLY_SATISFIABLE;
02068 else
02069 return false;
02070
02071 if (!(s >> str) || str != "tableau")
02072 return false;
02073
02074 if (!tableau.ascii_load(s))
02075 return false;
02076
02077 if (!(s >> str) || str != "working_cost(")
02078 return false;
02079
02080 dimension_type working_cost_dim;
02081
02082 if (!(s >> working_cost_dim))
02083 return false;
02084
02085 if (!(s >> str) || str != ")")
02086 return false;
02087
02088 if (!working_cost.ascii_load(s))
02089 return false;
02090
02091 if (!(s >> str) || str != "base(")
02092 return false;
02093
02094 dimension_type base_size;
02095 if (!(s >> base_size))
02096 return false;
02097
02098 if (!(s >> str) || str != ")")
02099 return false;
02100
02101 dimension_type base_value;
02102 for (dimension_type i = 0; i != base_size; ++i) {
02103 if (!(s >> base_value))
02104 return false;
02105 base.push_back(base_value);
02106 }
02107
02108 if (!(s >> str) || str != "last_generator")
02109 return false;
02110
02111 if (!last_generator.ascii_load(s))
02112 return false;
02113
02114 if (!(s >> str) || str != "mapping(")
02115 return false;
02116
02117 dimension_type mapping_size;
02118 if (!(s >> mapping_size))
02119 return false;
02120
02121 if (!(s >> str) || str != ")")
02122 return false;
02123
02124 dimension_type first_value;
02125 dimension_type second_value;
02126 dimension_type index;
02127
02128
02129
02130 if (tableau.num_columns() != 0)
02131 mapping.push_back(std::make_pair(0, 0));
02132
02133 for (dimension_type i = 1; i < mapping_size; ++i) {
02134 if (!(s >> index))
02135 return false;
02136 if (!(s >> str) || str != "->")
02137 return false;
02138 if (!(s >> first_value))
02139 return false;
02140 if (!(s >> str) || str != "->")
02141 return false;
02142 if (!(s >> second_value))
02143 return false;
02144 mapping.push_back(std::make_pair(first_value, second_value));
02145 }
02146
02147 if (!(s >> str) || str != "integer_variables")
02148 return false;
02149
02150 if (!i_variables.ascii_load(s))
02151 return false;
02152
02153 assert(OK());
02154 return true;
02155 }
02156
02158 std::ostream&
02159 PPL::IO_Operators::operator<<(std::ostream& s, const MIP_Problem& lp) {
02160 s << "Constraints:";
02161 for (MIP_Problem::const_iterator i = lp.constraints_begin(),
02162 i_end = lp.constraints_end(); i != i_end; ++i)
02163 s << "\n" << *i;
02164 s << "\nObjective function: "
02165 << lp.objective_function()
02166 << "\nOptimization mode: "
02167 << (lp.optimization_mode() == MAXIMIZATION
02168 ? "MAXIMIZATION"
02169 : "MINIMIZATION");
02170 s << "\nInteger variables: " << lp.integer_space_dimensions();
02171 return s;
02172 }