00001
00002
00003 #include <memory>
00004 #include <algorithm>
00005 #include <numeric>
00006 #include <cmath>
00007
00008 #include "CoinHelperFunctions.hpp"
00009 #include "CoinWarmStartBasis.hpp"
00010
00011 #include "BCP_problem_core.hpp"
00012 #include "BCP_lp_node.hpp"
00013 #include "BCP_lp_result.hpp"
00014 #include "BCP_lp_pool.hpp"
00015 #include "BCP_lp.hpp"
00016 #include "BCP_lp_branch.hpp"
00017 #include "BCP_lp_user.hpp"
00018 #include "BCP_lp_functions.hpp"
00019
00020
00021
00022
00023
00024
00025 inline static int BCP_compare_waiting_row_ptr(const BCP_lp_waiting_row* wrow0,
00026 const BCP_lp_waiting_row* wrow1){
00027 return wrow0->violation() < wrow1->violation();
00028 }
00029
00030
00031
00032
00033
00034
00035 inline static int BCP_compare_waiting_col_ptr(const BCP_lp_waiting_col* wcol0,
00036 const BCP_lp_waiting_col* wcol1){
00037 return wcol0->red_cost() > wcol1->red_cost();
00038 }
00039
00040
00041
00042 bool BCP_lp_fix_vars(BCP_lp_prob& p)
00043 {
00044 const BCP_lp_result& lpres = *p.lp_result;
00045
00046 BCP_var_set& vars = p.node->vars;
00047 const int varnum = vars.size();
00048 BCP_cut_set& cuts = p.node->cuts;
00049 const int cutnum = cuts.size();
00050
00051 int i;
00052 int newly_changed = 0;
00053
00054 BCP_lp_check_ub(p);
00055
00056 if (lpres.dj()) {
00057 p.user->reduced_cost_fixing(lpres.dj(), lpres.x(),
00058 p.ub() - lpres.objval() - p.granularity(),
00059 vars, newly_changed);
00060 }
00061 if (newly_changed > 0 && p.param(BCP_lp_par::LpVerb_VarTightening)) {
00062 printf("LP: Reduced cost fixing has changed the bounds on %i vars\n",
00063 newly_changed);
00064 }
00065
00066 p.var_bound_changes_since_logical_fixing += newly_changed;
00067
00068 BCP_vec<BCP_obj_status> var_status;
00069 var_status.reserve(varnum);
00070 for (i = 0; i < varnum; ++i)
00071 var_status.unchecked_push_back(vars[i]->status());
00072
00073 BCP_vec<BCP_obj_status> cut_status;
00074 cut_status.reserve(cutnum);
00075 for (i = 0; i < cutnum; ++i)
00076 cut_status.unchecked_push_back(cuts[i]->status());
00077
00078 bool lost_primal_feasibility = false;
00079
00080 BCP_vec<int> changed_pos;
00081 BCP_vec<double> new_bd;
00082
00083 p.user->logical_fixing(lpres, vars, cuts, var_status, cut_status,
00084 p.var_bound_changes_since_logical_fixing,
00085 changed_pos, new_bd);
00086
00087 if (2 * changed_pos.size() != new_bd.size())
00088 throw BCP_fatal_error("logical_fixing() returned uneven vectors\n");
00089
00090 const int change_num = changed_pos.size();
00091 if (change_num > 0){
00092 int bd_changes = 0;
00093 const double * x = lpres.x();
00094 const double petol = lpres.primalTolerance();
00095 for (i = 0; i < change_num; ++i) {
00096
00097
00098 const int pos = changed_pos[i];
00099 const double old_lb = vars[pos]->lb();
00100 const double old_ub = vars[pos]->ub();
00101 const double new_lb = new_bd[2*i];
00102 const double new_ub = new_bd[2*i+1];
00103 if (old_lb > new_lb + petol || old_ub < new_ub - petol) {
00104 throw BCP_fatal_error("logical fixing enlarged feas region!\n");
00105 }
00106
00107
00108 if (new_lb < old_lb + petol && new_ub > old_ub - petol) {
00109 printf("LP: BCP_lp_fix_vars(): WARNING: no change in bounds.\n");
00110 continue;
00111 }
00112
00113
00114
00115 new_bd[2*bd_changes] = new_lb;
00116 new_bd[2*bd_changes+1] = new_ub;
00117 changed_pos[bd_changes] = pos;
00118 ++bd_changes;
00119 if ((x[pos] < new_lb - petol) || (x[pos] > new_ub + petol))
00120 lost_primal_feasibility = true;
00121 }
00122 changed_pos.erase(changed_pos.entry(bd_changes), changed_pos.end());
00123 new_bd.erase(new_bd.entry(2*bd_changes), new_bd.end());
00124
00125 p.var_bound_changes_since_logical_fixing = 0;
00126
00127 if (bd_changes > 0) {
00128 p.lp_solver->setColSetBounds(changed_pos.begin(), changed_pos.end(),
00129 new_bd.begin());
00130 vars.set_lb_ub(changed_pos, new_bd.begin());
00131 if (bd_changes > 0 && p.param(BCP_lp_par::LpVerb_VarTightening)) {
00132 printf("LP: Logical fixing has changed the bounds on %i vars\n",
00133 bd_changes);
00134 }
00135 }
00136 }
00137
00138 return lost_primal_feasibility;
00139 }
00140
00141
00142
00143 static void
00144 BCP_lp_reset_positions(const BCP_vec<int>& deletable, BCP_vec<int>& pos,
00145 const bool error_if_deletable)
00146 {
00147 int i, j;
00148 const int delnum = deletable.size();
00149 const int posnum = pos.size();
00150 for (i = 0, j = 0; i < delnum && j < posnum; ) {
00151 #ifdef PARANOID
00152 if (error_if_deletable && deletable[i] == pos[j])
00153 throw BCP_fatal_error("Bad pos in BCP_lp_reset_positions()\n");
00154 #endif
00155 if (deletable[i] < pos[j]) {
00156 ++i;
00157 } else {
00158 pos[j] -= i;
00159 ++j;
00160 }
00161 }
00162 for ( ; j < posnum; ++j) {
00163 pos[j] -= delnum;
00164 }
00165 }
00166
00167
00168
00169 void
00170 BCP_delete_unwanted_candidates(const int num, const int added_num,
00171 const BCP_vec<int>* pos,
00172 BCP_vec<int>& deletable)
00173 {
00174 deletable.erase(std::lower_bound(deletable.begin(), deletable.end(),
00175 num-added_num), deletable.end());
00176 if (pos && !pos->empty()) {
00177 BCP_vec<int> ordered_pos(*pos);
00178 std::sort(ordered_pos.begin(), ordered_pos.end());
00179 const int size = ordered_pos.size();
00180 int i=num-added_num, j=0;
00181 while (i < num && j < size) {
00182 if (i < ordered_pos[j]) {
00183 deletable.push_back(i);
00184 ++i;
00185 continue;
00186 }
00187 if (i > ordered_pos[j]) {
00188 ++j;
00189 continue;
00190 }
00191 ++i;
00192 ++j;
00193 }
00194 for ( ; i < num; ++i) {
00195 deletable.push_back(i);
00196 }
00197 }
00198 }
00199
00200
00201 void BCP_lp_delete_cols_and_rows(BCP_lp_prob& p,
00202 BCP_lp_branching_object* can,
00203 const int added_colnum,
00204 const int added_rownum,
00205 const bool from_fathom,
00206 const bool force_delete)
00207 {
00208
00209
00210
00211
00212 int del_num;
00213
00214 BCP_lp_result& lpres = *p.lp_result;
00215 BCP_var_set& vars = p.node->vars;
00216 const size_t varnum = vars.size();
00217 BCP_cut_set& cuts = p.node->cuts;
00218 const size_t cutnum = cuts.size();
00219
00220 OsiSolverInterface* lp = p.lp_solver;
00221 CoinWarmStart* ws = lp->getWarmStart();
00222 CoinWarmStartBasis* bas = dynamic_cast<CoinWarmStartBasis*>(ws);
00223
00224 BCP_vec<int> deletable;
00225 deletable.reserve(CoinMax(varnum, cutnum));
00226
00227
00228 p.user->select_vars_to_delete(lpres, vars, cuts, from_fathom, deletable);
00229
00230
00231 del_num = deletable.size();
00232 if (del_num > 0) {
00233
00234 std::sort(deletable.begin(), deletable.end());
00235 deletable.erase(std::unique(deletable.begin(), deletable.end()),
00236 deletable.end());
00237 del_num = deletable.size();
00238 if (can && added_colnum) {
00239
00240 BCP_delete_unwanted_candidates(varnum, added_colnum,
00241 can->forced_var_pos, deletable);
00242 del_num = deletable.size();
00243 }
00244 if (bas) {
00245
00246
00247 BCP_vec<int> do_not_delete;
00248 for (int i = 0; i < del_num; ++i) {
00249 const int ind = deletable[i];
00250 if (can && ind >= (int)(varnum - added_colnum))
00251 continue;
00252 const CoinWarmStartBasis::Status stat = bas->getStructStatus(ind);
00253 if (stat == CoinWarmStartBasis::basic ||
00254 stat == CoinWarmStartBasis::isFree) {
00255
00256 do_not_delete.push_back(i);
00257 }
00258 }
00259 if (! do_not_delete.empty()) {
00260 deletable.unchecked_erase_by_index(do_not_delete);
00261 del_num = deletable.size();
00262 }
00263 }
00264 }
00265
00266 int min_by_frac =
00267 static_cast<int>(varnum * p.param(BCP_lp_par::DeletedColToCompress_Frac));
00268 int min_to_del = force_delete ?
00269 0 : CoinMax(p.param(BCP_lp_par::DeletedColToCompress_Min), min_by_frac);
00270
00271 if (del_num > min_to_del) {
00272 if (p.param(BCP_lp_par::LpVerb_MatrixCompression))
00273 printf("LP: Deleting %i columns from the matrix.\n", del_num);
00274 if (can) {
00275 if (can->forced_var_pos)
00276 BCP_lp_reset_positions(deletable, *can->forced_var_pos, true);
00277 if (can->implied_var_pos)
00278 BCP_lp_reset_positions(deletable, *can->implied_var_pos, false);
00279 }
00280 if (bas) {
00281 bas->deleteColumns(del_num, &deletable[0]);
00282 }
00283 p.lp_solver->deleteCols(del_num, &deletable[0]);
00284 purge_ptr_vector_by_index(dynamic_cast< BCP_vec<BCP_var*>& >(vars),
00285 deletable.begin(), deletable.end());
00286 p.local_cut_pool->rows_are_valid(false);
00287 }
00288
00289
00290 deletable.clear();
00291 p.user->select_cuts_to_delete(lpres, vars, cuts, from_fathom, deletable);
00292
00293
00294 del_num = deletable.size();
00295 if (del_num > 0) {
00296
00297 std::sort(deletable.begin(), deletable.end());
00298 deletable.erase(std::unique(deletable.begin(), deletable.end()),
00299 deletable.end());
00300 del_num = deletable.size();
00301 if (can && added_rownum) {
00302
00303 BCP_delete_unwanted_candidates(cutnum, added_rownum,
00304 can->forced_cut_pos, deletable);
00305 del_num = deletable.size();
00306 }
00307 if (bas) {
00308
00309
00310 BCP_vec<int> do_not_delete;
00311 for (int i = 0; i < del_num; ++i) {
00312 const int ind = deletable[i];
00313 if (can && ind >= (int)(cutnum - added_rownum))
00314 continue;
00315 const CoinWarmStartBasis::Status stat = bas->getArtifStatus(ind);
00316 if (stat != CoinWarmStartBasis::basic) {
00317
00318 do_not_delete.push_back(i);
00319 }
00320 }
00321 if (! do_not_delete.empty()) {
00322 deletable.unchecked_erase_by_index(do_not_delete);
00323 del_num = deletable.size();
00324 }
00325 }
00326 }
00327
00328 min_by_frac =
00329 static_cast<int>(cutnum * p.param(BCP_lp_par::DeletedRowToCompress_Frac));
00330 min_to_del = force_delete ?
00331 0 : CoinMax(p.param(BCP_lp_par::DeletedRowToCompress_Min), min_by_frac);
00332
00333 if (del_num > min_to_del) {
00334 if (p.param(BCP_lp_par::LpVerb_MatrixCompression))
00335 printf("LP: Deleting %i rows from the matrix.\n", del_num);
00336 if (can) {
00337 if (can->forced_cut_pos)
00338 BCP_lp_reset_positions(deletable, *can->forced_cut_pos, true);
00339 if (can->implied_cut_pos)
00340 BCP_lp_reset_positions(deletable, *can->implied_cut_pos, false);
00341 } else {
00342 if (p.param(BCP_lp_par::BranchOnCuts))
00343 p.node->cuts.move_deletable_to_pool(deletable, p.slack_pool);
00344 }
00345 if (bas) {
00346 bas->deleteRows(del_num, &deletable[0]);
00347 }
00348 p.lp_solver->deleteRows(deletable.size(), &deletable[0]);
00349 purge_ptr_vector_by_index(dynamic_cast< BCP_vec<BCP_cut*>& >(cuts),
00350 deletable.begin(), deletable.end());
00351 p.node->lb_at_cutgen.erase_by_index(deletable);
00352 p.local_var_pool->cols_are_valid(false);
00353 }
00354
00355 if (bas) {
00356 if (varnum != vars.size() || cutnum != cuts.size()) {
00357
00358 lp->setWarmStart(bas);
00359 }
00360 delete bas;
00361 }
00362 }
00363
00364
00365
00366 void BCP_lp_adjust_row_effectiveness(BCP_lp_prob& p)
00367 {
00368 if (p.param(BCP_lp_par::IneffectiveConstraints) != BCP_IneffConstr_None){
00369 const BCP_lp_result& lpres = *p.lp_result;
00370
00371 BCP_cut_set& cuts = p.node->cuts;
00372 int i, ineff = 0;
00373 const int cutnum = cuts.size();
00374
00375 switch (p.param(BCP_lp_par::IneffectiveConstraints)){
00376 case BCP_IneffConstr_None:
00377 break;
00378 case BCP_IneffConstr_NonzeroSlack:
00379 {
00380 const double * lhs = p.lp_result->lhs();
00381 const double * rlb = p.lp_solver->getRowLower();
00382 const double * rub = p.lp_solver->getRowUpper();
00383 const double petol = lpres.primalTolerance();
00384 for (i = p.core->cutnum(); i < cutnum; ++i) {
00385 BCP_cut *cut = cuts[i];
00386 if (rlb[i] + petol < lhs[i] && lhs[i] < rub[i] - petol) {
00387 cut->decrease_effective_count();
00388 if (! cut->is_non_removable())
00389 ++ineff;
00390 } else {
00391 cut->increase_effective_count();
00392 }
00393 }
00394 }
00395 break;
00396 case BCP_IneffConstr_ZeroDualValue:
00397 {
00398 const double * pi = p.lp_result->pi();
00399 const double detol = lpres.dualTolerance();
00400 for (i = p.core->cutnum(); i < cutnum; ++i) {
00401 BCP_cut *cut = cuts[i];
00402 if (CoinAbs(pi[i]) < detol) {
00403 cut->decrease_effective_count();
00404 if (! cut->is_non_removable())
00405 ++ineff;
00406 } else {
00407 cut->increase_effective_count();
00408 }
00409 }
00410 }
00411 break;
00412 }
00413 if (p.param(BCP_lp_par::LpVerb_RowEffectivenessCount))
00414 printf("LP: Row effectiveness: rownum: %i ineffective: %i\n",
00415 static_cast<int>(p.node->cuts.size()), ineff);
00416 }
00417 }
00418
00419
00420
00421 int BCP_lp_add_from_local_cut_pool(BCP_lp_prob& p)
00422 {
00423
00424 BCP_lp_cut_pool& cp = *p.local_cut_pool;
00425 size_t added_rows =
00426 std::min<size_t>(size_t(p.param(BCP_lp_par::MaxCutsAddedPerIteration)),
00427 cp.size());
00428
00429 if (added_rows == 0)
00430 return 0;
00431
00432 if (added_rows < cp.size()) {
00433
00434
00435
00436 switch (p.param(BCP_lp_par::CutViolationNorm)) {
00437 case BCP_CutViolationNorm_Plain:
00438 break;
00439 case BCP_CutViolationNorm_Distance:
00440 for (int i = cp.size() - 1; i >= 0; --i) {
00441 BCP_lp_waiting_row& cut = *cp[i];
00442 if (cut.violation() > 0) {
00443 cut.set_violation(cut.violation()/sqrt(cut.row()->normSquare()));
00444 }
00445 }
00446 break;
00447 case BCP_CutViolationNorm_Directional: {
00448 const double* c = p.lp_solver->getObjCoefficients();
00449 const int n = p.lp_solver->getNumCols();
00450 const double cnorm = sqrt(std::inner_product(c, c+n, c, 0.0));
00451 for (int i = cp.size() - 1; i >= 0; --i) {
00452 BCP_lp_waiting_row& cut = *cp[i];
00453 if (cut.violation() > 0) {
00454 cut.set_violation(cut.violation() * cnorm /
00455 fabs(cut.row()->dotProduct(c)));
00456 }
00457 }
00458 }
00459 break;
00460 default:
00461 throw BCP_fatal_error("LP: No violation norm specified!\n");
00462 }
00463
00464
00465 std::sort(cp.begin(), cp.end(), BCP_compare_waiting_row_ptr);
00466 }
00467
00468 #ifdef PARANOID
00469 if (! cp.rows_are_valid())
00470 throw BCP_fatal_error
00471 ("BCP_lp_add_from_local_cut_pool() : rows are not valid.\n");
00472 #endif
00473
00474
00475
00476 BCP_lp_cut_pool::const_iterator first = cp.entry(cp.size() - added_rows);
00477 BCP_lp_cut_pool::const_iterator last = cp.end();
00478 BCP_lp_cut_pool::const_iterator cpi = first - 1;
00479
00480
00481 const CoinPackedVectorBase** rows =
00482 new const CoinPackedVectorBase*[added_rows];
00483
00484 BCP_vec<double> rlb;
00485 rlb.reserve(added_rows);
00486
00487 BCP_vec<double> rub;
00488 rub.reserve(added_rows);
00489
00490 p.node->cuts.reserve(p.node->cuts.size() + added_rows);
00491
00492 cpi = first - 1;
00493 while (++cpi != last) {
00494 BCP_row * row = (*cpi)->row();
00495 BCP_cut * cut = (*cpi)->cut();
00496 *rows++ = row;
00497 rlb.unchecked_push_back(row->LowerBound());
00498 rub.unchecked_push_back(row->UpperBound());
00499 cut->set_effective_count(0);
00500 p.node->cuts.unchecked_push_back(cut);
00501 }
00502 rows -= added_rows;
00503 p.lp_solver->addRows(added_rows, rows, rlb.begin(), rub.begin());
00504 cpi = first - 1;
00505 while (++cpi != last) {
00506 (*cpi)->clear_cut();
00507 }
00508
00509 delete[] rows;
00510
00511 int leftover = cp.size() - added_rows;
00512 const double frac = p.param(BCP_lp_par::MaxLeftoverCutFrac);
00513 const int maxnum = p.param(BCP_lp_par::MaxLeftoverCutNum);
00514 if (frac < 1.0) {
00515 leftover = (frac <= 0.0) ? 0 : static_cast<int>(frac * leftover);
00516 }
00517 if (leftover > maxnum)
00518 leftover = maxnum > 0 ? maxnum : 0;
00519
00520 purge_ptr_vector(dynamic_cast<BCP_vec<BCP_lp_waiting_row*>&>(cp),
00521 cp.entry(leftover), cp.end());
00522
00523 p.node->lb_at_cutgen.insert(p.node->lb_at_cutgen.end(), added_rows,
00524 p.lp_result->objval());
00525 return added_rows;
00526 }
00527
00528
00529
00530 int BCP_lp_add_from_local_var_pool(BCP_lp_prob& p)
00531 {
00532
00533 BCP_lp_var_pool& vp = *p.local_var_pool;
00534 size_t added_cols =
00535 std::min<size_t>(size_t(p.param(BCP_lp_par::MaxVarsAddedPerIteration)),
00536 vp.size());
00537
00538 if (added_cols == 0)
00539 return 0;
00540
00541 if (added_cols < vp.size()) {
00542
00543
00544 std::sort(vp.begin(), vp.end(), BCP_compare_waiting_col_ptr);
00545 }
00546
00547 #ifdef PARANOID
00548 if (! vp.cols_are_valid())
00549 throw BCP_fatal_error
00550 ("BCP_lp_add_from_local_var_pool() : cols are not valid.\n");
00551 #endif
00552
00553
00554
00555 BCP_lp_var_pool::const_iterator first = vp.entry(vp.size() - added_cols);
00556 BCP_lp_var_pool::const_iterator last = vp.end();
00557 BCP_lp_var_pool::const_iterator vpi = first - 1;
00558
00559
00560 const CoinPackedVectorBase** cols =
00561 new const CoinPackedVectorBase*[added_cols];
00562
00563 BCP_vec<double> clb;
00564 clb.reserve(added_cols);
00565
00566 BCP_vec<double> cub;
00567 cub.reserve(added_cols);
00568
00569 BCP_vec<double> obj;
00570 obj.reserve(added_cols);
00571
00572 p.node->vars.reserve(p.node->vars.size() + added_cols);
00573
00574 vpi = first - 1;
00575 while (++vpi != last) {
00576 BCP_col* col = (*vpi)->col();
00577 BCP_var* var = (*vpi)->var();
00578 p.node->vars.unchecked_push_back(var);
00579
00580 if (var->var_type() != BCP_ContinuousVar) {
00581 var->set_lb(floor(var->lb()+1e-8));
00582 var->set_ub(ceil(var->ub()-1e-8));
00583 }
00584 *cols++ = col;
00585 clb.unchecked_push_back(col->LowerBound());
00586 cub.unchecked_push_back(col->UpperBound());
00587 obj.unchecked_push_back(col->Objective());
00588 }
00589 cols -= added_cols;
00590 p.lp_solver->addCols(added_cols, cols,
00591 clb.begin(), cub.begin(), obj.begin());
00592 vpi = first - 1;
00593 while (++vpi != last) {
00594 (*vpi)->clear_var();
00595 }
00596
00597 delete[] cols;
00598
00599 purge_ptr_vector(dynamic_cast<BCP_vec<BCP_lp_waiting_col*>&>(vp),
00600 vp.entry(vp.size() - added_cols), vp.end());
00601
00602 return added_cols;
00603 }
00604
00605