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 for (int i = num - added_num; i < num; ++i) {
00177 deletable.push_back(i);
00178 }
00179 if (pos && !pos->empty()) {
00180 BCP_vec<int>::iterator ifirst =
00181 std::lower_bound(deletable.begin(), deletable.end(), (*pos)[0]);
00182 BCP_vec<int>::iterator ilast =
00183 std::upper_bound(ifirst, deletable.end(), pos->back());
00184 deletable.erase(ifirst, ilast);
00185 }
00186 }
00187
00188
00189 void BCP_lp_delete_cols_and_rows(BCP_lp_prob& p,
00190 BCP_lp_branching_object* can,
00191 const int added_colnum,
00192 const int added_rownum,
00193 const bool from_fathom,
00194 const bool force_delete)
00195 {
00196
00197
00198
00199
00200 int del_num;
00201
00202 BCP_lp_result& lpres = *p.lp_result;
00203 BCP_var_set& vars = p.node->vars;
00204 const size_t varnum = vars.size();
00205 BCP_cut_set& cuts = p.node->cuts;
00206 const size_t cutnum = cuts.size();
00207
00208 OsiSolverInterface* lp = p.lp_solver;
00209 CoinWarmStart* ws = lp->getWarmStart();
00210 CoinWarmStartBasis* bas = dynamic_cast<CoinWarmStartBasis*>(ws);
00211
00212 BCP_vec<int> deletable;
00213 deletable.reserve(CoinMax(varnum, cutnum));
00214
00215
00216 p.user->select_vars_to_delete(lpres, vars, cuts, from_fathom, deletable);
00217
00218
00219 del_num = deletable.size();
00220 if (del_num > 0) {
00221
00222 std::sort(deletable.begin(), deletable.end());
00223 deletable.erase(std::unique(deletable.begin(), deletable.end()),
00224 deletable.end());
00225 del_num = deletable.size();
00226 if (can && added_colnum) {
00227
00228 BCP_delete_unwanted_candidates(varnum, added_colnum,
00229 can->forced_var_pos, deletable);
00230 del_num = deletable.size();
00231 }
00232 if (bas) {
00233
00234
00235 BCP_vec<int> do_not_delete;
00236 for (int i = 0; i < del_num; ++i) {
00237 const int ind = deletable[i];
00238 if (can && ind >= (int)(varnum - added_colnum))
00239 continue;
00240 const CoinWarmStartBasis::Status stat = bas->getStructStatus(ind);
00241 if (stat == CoinWarmStartBasis::basic ||
00242 stat == CoinWarmStartBasis::isFree) {
00243
00244 do_not_delete.push_back(i);
00245 }
00246 }
00247 if (! do_not_delete.empty()) {
00248 deletable.unchecked_erase_by_index(do_not_delete);
00249 del_num = deletable.size();
00250 }
00251 }
00252 }
00253
00254 int min_by_frac =
00255 static_cast<int>(varnum * p.param(BCP_lp_par::DeletedColToCompress_Frac));
00256 int min_to_del = force_delete ?
00257 0 : CoinMax(p.param(BCP_lp_par::DeletedColToCompress_Min), min_by_frac);
00258
00259 if (del_num > min_to_del) {
00260 if (p.param(BCP_lp_par::LpVerb_MatrixCompression))
00261 printf("LP: Deleting %i columns from the matrix.\n", del_num);
00262 if (can) {
00263 if (can->forced_var_pos)
00264 BCP_lp_reset_positions(deletable, *can->forced_var_pos, true);
00265 if (can->implied_var_pos)
00266 BCP_lp_reset_positions(deletable, *can->implied_var_pos, false);
00267 }
00268 if (bas) {
00269 bas->deleteColumns(del_num, &deletable[0]);
00270 }
00271 p.lp_solver->deleteCols(del_num, &deletable[0]);
00272 purge_ptr_vector_by_index(dynamic_cast< BCP_vec<BCP_var*>& >(vars),
00273 deletable.begin(), deletable.end());
00274 p.local_cut_pool->rows_are_valid(false);
00275 }
00276
00277
00278 deletable.clear();
00279 p.user->select_cuts_to_delete(lpres, vars, cuts, from_fathom, deletable);
00280
00281
00282 del_num = deletable.size();
00283 if (del_num > 0) {
00284
00285 std::sort(deletable.begin(), deletable.end());
00286 deletable.erase(std::unique(deletable.begin(), deletable.end()),
00287 deletable.end());
00288 del_num = deletable.size();
00289 if (can && added_rownum) {
00290
00291 BCP_delete_unwanted_candidates(cutnum, added_rownum,
00292 can->forced_cut_pos, deletable);
00293 del_num = deletable.size();
00294 }
00295 if (bas) {
00296
00297
00298 BCP_vec<int> do_not_delete;
00299 for (int i = 0; i < del_num; ++i) {
00300 const int ind = deletable[i];
00301 if (can && ind >= (int)(cutnum - added_rownum))
00302 continue;
00303 const CoinWarmStartBasis::Status stat = bas->getArtifStatus(ind);
00304 if (stat != CoinWarmStartBasis::basic) {
00305
00306 do_not_delete.push_back(i);
00307 }
00308 }
00309 if (! do_not_delete.empty()) {
00310 deletable.unchecked_erase_by_index(do_not_delete);
00311 del_num = deletable.size();
00312 }
00313 }
00314 }
00315
00316 min_by_frac =
00317 static_cast<int>(cutnum * p.param(BCP_lp_par::DeletedRowToCompress_Frac));
00318 min_to_del = force_delete ?
00319 0 : CoinMax(p.param(BCP_lp_par::DeletedRowToCompress_Min), min_by_frac);
00320
00321 if (del_num > min_to_del) {
00322 if (p.param(BCP_lp_par::LpVerb_MatrixCompression))
00323 printf("LP: Deleting %i rows from the matrix.\n", del_num);
00324 if (can) {
00325 if (can->forced_cut_pos)
00326 BCP_lp_reset_positions(deletable, *can->forced_cut_pos, true);
00327 if (can->implied_cut_pos)
00328 BCP_lp_reset_positions(deletable, *can->implied_cut_pos, false);
00329 } else {
00330 if (p.param(BCP_lp_par::BranchOnCuts))
00331 p.node->cuts.move_deletable_to_pool(deletable, p.slack_pool);
00332 }
00333 if (bas) {
00334 bas->deleteRows(del_num, &deletable[0]);
00335 }
00336 p.lp_solver->deleteRows(deletable.size(), &deletable[0]);
00337 purge_ptr_vector_by_index(dynamic_cast< BCP_vec<BCP_cut*>& >(cuts),
00338 deletable.begin(), deletable.end());
00339 p.node->lb_at_cutgen.erase_by_index(deletable);
00340 p.local_var_pool->cols_are_valid(false);
00341 }
00342
00343 if (bas) {
00344 if (varnum != vars.size() || cutnum != cuts.size()) {
00345
00346 lp->setWarmStart(bas);
00347 }
00348 delete bas;
00349 }
00350 }
00351
00352
00353
00354 void BCP_lp_adjust_row_effectiveness(BCP_lp_prob& p)
00355 {
00356 if (p.param(BCP_lp_par::IneffectiveConstraints) != BCP_IneffConstr_None){
00357 const BCP_lp_result& lpres = *p.lp_result;
00358
00359 BCP_cut_set& cuts = p.node->cuts;
00360 int i, ineff = 0;
00361 const int cutnum = cuts.size();
00362
00363 switch (p.param(BCP_lp_par::IneffectiveConstraints)){
00364 case BCP_IneffConstr_None:
00365 break;
00366 case BCP_IneffConstr_NonzeroSlack:
00367 {
00368 const double * lhs = p.lp_result->lhs();
00369 const double * rlb = p.lp_solver->getRowLower();
00370 const double * rub = p.lp_solver->getRowUpper();
00371 const double petol = lpres.primalTolerance();
00372 for (i = p.core->cutnum(); i < cutnum; ++i) {
00373 BCP_cut *cut = cuts[i];
00374 if (rlb[i] + petol < lhs[i] && lhs[i] < rub[i] - petol) {
00375 cut->decrease_effective_count();
00376 if (! cut->is_non_removable())
00377 ++ineff;
00378 } else {
00379 cut->increase_effective_count();
00380 }
00381 }
00382 }
00383 break;
00384 case BCP_IneffConstr_ZeroDualValue:
00385 {
00386 const double * pi = p.lp_result->pi();
00387 const double detol = lpres.dualTolerance();
00388 for (i = p.core->cutnum(); i < cutnum; ++i) {
00389 BCP_cut *cut = cuts[i];
00390 if (CoinAbs(pi[i]) < detol) {
00391 cut->decrease_effective_count();
00392 if (! cut->is_non_removable())
00393 ++ineff;
00394 } else {
00395 cut->increase_effective_count();
00396 }
00397 }
00398 }
00399 break;
00400 }
00401 if (p.param(BCP_lp_par::LpVerb_RowEffectivenessCount))
00402 printf("LP: Row effectiveness: rownum: %i ineffective: %i\n",
00403 static_cast<int>(p.node->cuts.size()), ineff);
00404 }
00405 }
00406
00407
00408
00409 int BCP_lp_add_from_local_cut_pool(BCP_lp_prob& p)
00410 {
00411
00412 BCP_lp_cut_pool& cp = *p.local_cut_pool;
00413 size_t added_rows =
00414 std::min<size_t>(size_t(p.param(BCP_lp_par::MaxCutsAddedPerIteration)),
00415 cp.size());
00416
00417 if (added_rows == 0)
00418 return 0;
00419
00420 if (added_rows < cp.size()) {
00421
00422
00423
00424 switch (p.param(BCP_lp_par::CutViolationNorm)) {
00425 case BCP_CutViolationNorm_Plain:
00426 break;
00427 case BCP_CutViolationNorm_Distance:
00428 for (int i = cp.size() - 1; i >= 0; --i) {
00429 BCP_lp_waiting_row& cut = *cp[i];
00430 if (cut.violation() > 0) {
00431 cut.set_violation(cut.violation()/sqrt(cut.row()->normSquare()));
00432 }
00433 }
00434 break;
00435 case BCP_CutViolationNorm_Directional: {
00436 const double* c = p.lp_solver->getObjCoefficients();
00437 const int n = p.lp_solver->getNumCols();
00438 const double cnorm = sqrt(std::inner_product(c, c+n, c, 0.0));
00439 for (int i = cp.size() - 1; i >= 0; --i) {
00440 BCP_lp_waiting_row& cut = *cp[i];
00441 if (cut.violation() > 0) {
00442 cut.set_violation(cut.violation() * cnorm /
00443 fabs(cut.row()->dotProduct(c)));
00444 }
00445 }
00446 }
00447 break;
00448 default:
00449 throw BCP_fatal_error("LP: No violation norm specified!\n");
00450 }
00451
00452
00453 std::sort(cp.begin(), cp.end(), BCP_compare_waiting_row_ptr);
00454 }
00455
00456 #ifdef PARANOID
00457 if (! cp.rows_are_valid())
00458 throw BCP_fatal_error
00459 ("BCP_lp_add_from_local_cut_pool() : rows are not valid.\n");
00460 #endif
00461
00462
00463
00464 BCP_lp_cut_pool::const_iterator first = cp.entry(cp.size() - added_rows);
00465 BCP_lp_cut_pool::const_iterator last = cp.end();
00466 BCP_lp_cut_pool::const_iterator cpi = first - 1;
00467
00468
00469 const CoinPackedVectorBase** rows =
00470 new const CoinPackedVectorBase*[added_rows];
00471
00472 BCP_vec<double> rlb;
00473 rlb.reserve(added_rows);
00474
00475 BCP_vec<double> rub;
00476 rub.reserve(added_rows);
00477
00478 p.node->cuts.reserve(p.node->cuts.size() + added_rows);
00479
00480 cpi = first - 1;
00481 while (++cpi != last) {
00482 BCP_row * row = (*cpi)->row();
00483 BCP_cut * cut = (*cpi)->cut();
00484 *rows++ = row;
00485 rlb.unchecked_push_back(row->LowerBound());
00486 rub.unchecked_push_back(row->UpperBound());
00487 cut->set_effective_count(0);
00488 p.node->cuts.unchecked_push_back(cut);
00489 }
00490 rows -= added_rows;
00491 p.lp_solver->addRows(added_rows, rows, rlb.begin(), rub.begin());
00492 cpi = first - 1;
00493 while (++cpi != last) {
00494 (*cpi)->clear_cut();
00495 }
00496
00497 delete[] rows;
00498
00499 int leftover = cp.size() - added_rows;
00500 const double frac = p.param(BCP_lp_par::MaxLeftoverCutFrac);
00501 const int maxnum = p.param(BCP_lp_par::MaxLeftoverCutNum);
00502 if (frac < 1.0) {
00503 leftover = (frac <= 0.0) ? 0 : static_cast<int>(frac * leftover);
00504 }
00505 if (leftover > maxnum)
00506 leftover = maxnum > 0 ? maxnum : 0;
00507
00508 purge_ptr_vector(dynamic_cast<BCP_vec<BCP_lp_waiting_row*>&>(cp),
00509 cp.entry(leftover), cp.end());
00510
00511 p.node->lb_at_cutgen.insert(p.node->lb_at_cutgen.end(), added_rows,
00512 p.lp_result->objval());
00513 return added_rows;
00514 }
00515
00516
00517
00518 int BCP_lp_add_from_local_var_pool(BCP_lp_prob& p)
00519 {
00520
00521 BCP_lp_var_pool& vp = *p.local_var_pool;
00522 size_t added_cols =
00523 std::min<size_t>(size_t(p.param(BCP_lp_par::MaxVarsAddedPerIteration)),
00524 vp.size());
00525
00526 if (added_cols == 0)
00527 return 0;
00528
00529 if (added_cols < vp.size()) {
00530
00531
00532 std::sort(vp.begin(), vp.end(), BCP_compare_waiting_col_ptr);
00533 }
00534
00535 #ifdef PARANOID
00536 if (! vp.cols_are_valid())
00537 throw BCP_fatal_error
00538 ("BCP_lp_add_from_local_var_pool() : cols are not valid.\n");
00539 #endif
00540
00541
00542
00543 BCP_lp_var_pool::const_iterator first = vp.entry(vp.size() - added_cols);
00544 BCP_lp_var_pool::const_iterator last = vp.end();
00545 BCP_lp_var_pool::const_iterator vpi = first - 1;
00546
00547
00548 const CoinPackedVectorBase** cols =
00549 new const CoinPackedVectorBase*[added_cols];
00550
00551 BCP_vec<double> clb;
00552 clb.reserve(added_cols);
00553
00554 BCP_vec<double> cub;
00555 cub.reserve(added_cols);
00556
00557 BCP_vec<double> obj;
00558 obj.reserve(added_cols);
00559
00560 p.node->vars.reserve(p.node->vars.size() + added_cols);
00561
00562 vpi = first - 1;
00563 while (++vpi != last) {
00564 BCP_col* col = (*vpi)->col();
00565 BCP_var* var = (*vpi)->var();
00566 p.node->vars.unchecked_push_back(var);
00567
00568 if (var->var_type() != BCP_ContinuousVar) {
00569 var->set_lb(floor(var->lb()+1e-8));
00570 var->set_ub(ceil(var->ub()-1e-8));
00571 }
00572 *cols++ = col;
00573 clb.unchecked_push_back(col->LowerBound());
00574 cub.unchecked_push_back(col->UpperBound());
00575 obj.unchecked_push_back(col->Objective());
00576 }
00577 cols -= added_cols;
00578 p.lp_solver->addCols(added_cols, cols,
00579 clb.begin(), cub.begin(), obj.begin());
00580 vpi = first - 1;
00581 while (++vpi != last) {
00582 (*vpi)->clear_var();
00583 }
00584
00585 delete[] cols;
00586
00587 purge_ptr_vector(dynamic_cast<BCP_vec<BCP_lp_waiting_col*>&>(vp),
00588 vp.entry(vp.size() - added_cols), vp.end());
00589
00590 return added_cols;
00591 }
00592
00593