00001
00002
00003
00004 #include <cstdio>
00005 #include <algorithm>
00006
00007 #include "CoinTime.hpp"
00008
00009 #include "OsiVolSolverInterface.hpp"
00010 #include "OsiClpSolverInterface.hpp"
00011 #include "ClpDualRowSteepest.hpp"
00012 #include "ClpPrimalColumnSteepest.hpp"
00013 #include "ClpSimplex.hpp"
00014
00015 #include "MC_lp.hpp"
00016
00017 #include "BCP_math.hpp"
00018 #include "BCP_lp.hpp"
00019 #include "BCP_lp_node.hpp"
00020 #include "BCP_lp_functions.hpp"
00021
00022
00023
00024 bool
00025 MC_lp::is_gap_tailoff_rel(const int k, const double minimp,
00026 const double objval) const
00027 {
00028 const int iteration_count = current_iteration();
00029 if (k <= 0 || iteration_count <= k)
00030 return false;
00031
00032 const int hist_size = CoinMin(iteration_count, hist_len);
00033 const double* obj_k_hist = objhist + (hist_size - k);
00034
00035
00036
00037 const double ub = upper_bound();
00038 return ((ub - objval) > (ub - obj_k_hist[0]) * (1-minimp));
00039 }
00040
00041
00042
00043 bool
00044 MC_lp::is_lb_tailoff_rel(const int k, const double minimp,
00045 const double objval) const
00046 {
00047 const int iteration_count = current_iteration();
00048 if (k <= 0 || iteration_count <= k)
00049 return false;
00050
00051 const int hist_size = CoinMin(iteration_count, hist_len);
00052 const double* obj_k_hist = objhist + (hist_size - k);
00053
00054
00055 return (objval - obj_k_hist[0] < minimp * CoinAbs(obj_k_hist[0]));
00056 }
00057
00058
00059
00060 bool
00061 MC_lp::is_lb_tailoff_abs(const int k, const double minimp,
00062 const double objval) const
00063 {
00064
00065
00066 const int iteration_count = current_iteration();
00067 if (k <= 0 || iteration_count <= k)
00068 return false;
00069
00070 const int hist_size = CoinMin(iteration_count, hist_len);
00071 const double* obj_k_hist = objhist + (hist_size - k);
00072
00073
00074 return (objval - obj_k_hist[0] < minimp);
00075 }
00076
00077
00078
00079 void
00080 MC_lp::tailoff_test(bool& tailoff_gap_rel, bool& tailoff_lb_abs,
00081 bool& tailoff_lb_rel, const double objval) const
00082 {
00083 tailoff_gap_rel =
00084 is_gap_tailoff_rel(par.entry(MC_lp_par::TailoffGapRelMinItcount),
00085 par.entry(MC_lp_par::TailoffGapRelMinImprovement),
00086 objval);
00087 tailoff_lb_abs =
00088 is_lb_tailoff_abs(par.entry(MC_lp_par::TailoffLbAbsMinItcount),
00089 par.entry(MC_lp_par::TailoffLbAbsMinImprovement),
00090 objval);
00091 tailoff_lb_rel =
00092 is_lb_tailoff_rel(par.entry(MC_lp_par::TailoffLbRelMinItcount),
00093 par.entry(MC_lp_par::TailoffLbRelMinImprovement),
00094 objval);
00095 const int iteration_count = current_iteration();
00096 const int hist_size = CoinMin(iteration_count - 1, hist_len);
00097 printf("MC: Tailoff check: objval: %.3f, UB: %.0f\n",
00098 objval, upper_bound());
00099 for (int i = hist_size - 1; i >= 0; --i) {
00100 printf(" LB[%3i]: %.3f\n", i-hist_size, objhist[i]);
00101 }
00102 printf(" gap_rel: %i lb_abs: %i lb_rel: %i\n",
00103 tailoff_gap_rel ? 1 : 0,
00104 tailoff_lb_abs ? 1 : 0,
00105 tailoff_lb_rel ? 1 : 0);
00106 }
00107
00108
00109
00110
00111 bool
00112 MC_solveClp(ClpSimplex& clp_simplex, OsiVolSolverInterface* vollp)
00113 {
00114 double t = CoinCpuTime();
00115
00116
00117
00118
00119 clp_simplex.scaling(0);
00120 ClpDualRowSteepest steep;
00121 clp_simplex.setDualRowPivotAlgorithm(steep);
00122 ClpPrimalColumnSteepest steepP;
00123 clp_simplex.setPrimalColumnPivotAlgorithm(steepP);
00124
00125 ClpSolve clp_options;
00126
00127 clp_options.setSolveType(ClpSolve::useDual);
00128
00129 clp_options.setPresolveType(ClpSolve::presolveOn);
00150 clp_options.setSpecialOption(0, 1);
00151 clp_options.setSpecialOption(2, 0);
00152 clp_options.setSpecialOption(3, 0);
00153 clp_simplex.initialSolve(clp_options);
00154
00155 t = CoinCpuTime() - t;
00156 printf("LP: exact optimization took %.3f seconds\n", t);
00157
00158 if (clp_simplex.isProvenPrimalInfeasible()) {
00159 printf("MC: Solving LP to optimality: infeas.\n");
00160 printf("MC: Forcing BCP to prune the node.\n");
00161 return false;
00162 }
00163 return true;
00164 }
00165
00166
00167
00168 OsiSolverInterface*
00169 MC_lp::solveToOpt(OsiVolSolverInterface* vollp,
00170 const BCP_lp_result& lpres,
00171 const BCP_vec<BCP_var*>& vars,
00172 const BCP_vec<BCP_cut*>& cuts,
00173 double& exact_obj)
00174 {
00175 if ((par.entry(MC_lp_par::LpSolver) & MC_UseClp) == 0) {
00176 throw BCP_fatal_error("\
00177 MC_lp::solveToOpt: got here but no simplex based solver is enabled in the\n\
00178 MC_LpSolver parameter.\n");
00179 }
00180
00181 OsiClpSolverInterface* solver = new OsiClpSolverInterface;
00182 solver->loadProblem(*vollp->getMatrixByCol(), vollp->getColLower(),
00183 vollp->getColUpper(), vollp->getObjCoefficients(),
00184 vollp->getRowLower(), vollp->getRowUpper());
00185 ClpSimplex& clp = *solver->getModelPtr();
00186
00187 if (par.entry(MC_lp_par::ExplicitSlacksInOpt)) {
00188 ClpSimplex clp_slack(clp);
00189 ClpSimplex& model = clp;
00190 ClpSimplex& model2 = clp_slack;
00191
00192 int numcols = model.numberColumns();
00193 int numrows = model.numberRows();
00194
00195 int * addStarts = new int [numrows+1];
00196 int * addRow = new int[numrows];
00197 double * addElement = new double[numrows];
00198 double * newUpper = new double[numrows];
00199 double * newLower = new double[numrows];
00200
00201 double * lower = model2.rowLower();
00202 double * upper = model2.rowUpper();
00203 int iRow, iColumn;
00204
00205
00206 for (iRow=0;iRow<numrows;iRow++) {
00207 newUpper[iRow]=upper[iRow];
00208 upper[iRow]=0.0;
00209 newLower[iRow]=lower[iRow];
00210 lower[iRow]=0.0;
00211 addRow[iRow]=iRow;
00212 addElement[iRow]=-1.0;
00213 addStarts[iRow]=iRow;
00214 }
00215 addStarts[numrows]=numrows;
00216 model2.addColumns(numrows,newLower,newUpper,NULL,
00217 addStarts,addRow,addElement);
00218 delete [] addStarts;
00219 delete [] addRow;
00220 delete [] addElement;
00221 delete [] newLower;
00222 delete [] newUpper;
00223
00224
00225 model2.transposeTimes(1.0,vollp->getRowPrice(),model2.objective());
00226
00227 if (!MC_solveClp(model2, vollp)) {
00228 delete solver;
00229 return NULL;
00230 }
00231 int lookup[]={0,1,2,3,0,3};
00232 CoinWarmStartBasis basis;
00233 basis.setSize(numcols,numrows);
00234 for (iColumn = 0; iColumn < numcols; iColumn++) {
00235 int iStatus = model2.getColumnStatus(iColumn);
00236 iStatus = lookup[iStatus];
00237 basis.setStructStatus(iColumn,(CoinWarmStartBasis::Status) iStatus);
00238 }
00239 for (iRow=0;iRow<numrows;iRow++) {
00240 int iStatus = model2.getRowStatus(iRow) == ClpSimplex::basic ?
00241 ClpSimplex::basic : model2.getColumnStatus(iRow+numcols);
00242 iStatus = lookup[iStatus];
00243 basis.setArtifStatus(iRow,(CoinWarmStartBasis::Status) iStatus);
00244 }
00245 solver->setWarmStart(&basis);
00246 solver->resolve();
00247 } else {
00248 if (!MC_solveClp(clp, vollp)) {
00249 delete solver;
00250 return NULL;
00251 }
00252 }
00253
00254 exact_obj = solver->getObjValue();
00255
00256 soln = mc_generate_heuristic_solution(solver->getColSolution(), vars, cuts);
00257 if (soln && soln->objective_value() < upper_bound()) {
00258 send_feasible_solution(soln);
00259 }
00260
00261
00262 printf("MC: Solving LP to optimality before branching: %.3f",
00263 exact_obj);
00264 const double gap =
00265 upper_bound() - exact_obj - get_param(BCP_lp_par::Granularity);
00266 printf(" Gap: %.3f\n",gap);
00267 if (gap < 0) {
00268 printf("MC: exact LP solving proved fathomability.\n");
00269 delete solver;
00270 return NULL;
00271 }
00272
00273
00274 int fix = 0;
00275 reduced_cost_fixing(solver->getReducedCost(), solver->getColSolution(), gap,
00276
00277 getLpProblemPointer()->node->vars, fix);
00278 printf("MC: rc fixing with exact opt and dj of exact fixed %i vars.\n", fix);
00279 return solver;
00280 }
00281
00282
00283
00284 void
00285 MC_lp::choose_branching_vars(const BCP_vec<BCP_var*>& vars, const double * x,
00286 const int cand_num,
00287 BCP_vec<BCP_lp_branching_object*>& cands)
00288 {
00289
00290
00291 const int m = mc.num_edges;
00292 const MC_graph_edge* edges = mc.edges;
00293
00294 BCP_vec<int> perm;
00295 perm.reserve(m);
00296
00297 BCP_vec<double> dist;
00298 dist.reserve(m);
00299
00300 const double etol = par.entry(MC_lp_par::IntegerTolerance);
00301
00302 int i;
00303
00304
00305
00306 if (mc.ising_triangles) {
00307 const int grid = static_cast<int>(sqrt(mc.num_nodes + 1.0));
00308 i = 2 * grid * grid;
00309 } else {
00310 i = 0;
00311 }
00312 for ( ; i < m; ++i) {
00313 const double xi = x[i];
00314 if (xi > etol && xi < 1-etol) {
00315 perm.unchecked_push_back(i);
00316 const double w = (0.5 - CoinAbs(xi-0.5)) * edges[i].cost;
00317 dist.unchecked_push_back(-CoinAbs(w));
00318 }
00319 }
00320
00321 const int distsize = dist.size();
00322 if (distsize > cand_num) {
00323 CoinSort_2(dist.begin(), dist.end(), perm.begin());
00324 perm.erase(perm.entry(cand_num), perm.end());
00325 }
00326
00327 if (perm.size() > 0) {
00328 append_branching_vars(x, vars, perm, cands);
00329 }
00330 }
00331
00332
00333
00334 BCP_branching_decision
00335 MC_lp::select_branching_candidates(const BCP_lp_result& lpres,
00336 const BCP_vec<BCP_var*>& vars,
00337 const BCP_vec<BCP_cut*>& cuts,
00338 const BCP_lp_var_pool& local_var_pool,
00339 const BCP_lp_cut_pool& local_cut_pool,
00340 BCP_vec<BCP_lp_branching_object*>& cands,
00341 bool force_branch)
00342 {
00343 OsiSolverInterface* lp = getLpProblemPointer()->lp_solver;
00344 OsiVolSolverInterface* vollp = dynamic_cast<OsiVolSolverInterface*>(lp);
00345 OsiSolverInterface* exact_solver = NULL;
00346 double objval = lpres.objval();
00347 bool tailoff_gap_rel = false;
00348 bool tailoff_lb_abs = false;
00349 bool tailoff_lb_rel = false;
00350
00351 const int solver_t = par.entry(MC_lp_par::LpSolver);
00352 const bool has_vol = (solver_t & MC_UseVol);
00353 const bool has_simplex = (solver_t & MC_UseClp);
00354
00355 int pool_size = local_cut_pool.size();
00356
00357 if (current_index() != 0 && has_vol && has_simplex) {
00358 started_exact = true;
00359 }
00360
00361 const int iteration_count = current_iteration();
00362 bool solve_to_opt = false;
00363
00364 if (started_exact) {
00365
00366 solve_to_opt = true;
00367 } else {
00368 tailoff_test(tailoff_gap_rel, tailoff_lb_abs, tailoff_lb_rel, objval);
00369 if (has_vol && has_simplex &&
00370 (tailoff_gap_rel && (tailoff_lb_abs || tailoff_lb_rel))) {
00371 solve_to_opt = true;
00372 }
00373 }
00374
00375 if (solve_to_opt) {
00376 exact_solver = solveToOpt(vollp, lpres, vars, cuts, objval);
00377 if (par.entry(MC_lp_par::OnceOptAlwaysOpt)) {
00378 started_exact = true;
00379 }
00380 if (exact_solver == NULL)
00381 return BCP_DoNotBranch_Fathomed;
00382 tailoff_test(tailoff_gap_rel, tailoff_lb_abs, tailoff_lb_rel, objval);
00383
00384
00385 BCP_lp_prob* p = getLpProblemPointer();
00386 BCP_lp_cut_pool& cp = *p->local_cut_pool;
00387 BCP_vec<BCP_cut*> new_cuts;
00388 BCP_vec<BCP_row*> new_rows;
00389 generate_cuts_in_lp(exact_solver->getColSolution(),
00390 exact_solver->getRowActivity(),
00391 exact_solver->getObjValue(),
00392 vars, cuts, new_cuts, new_rows);
00393 const int new_size = new_cuts.size();
00394 if (new_rows.size() != 0) {
00395 cp.reserve(cp.size() + new_size);
00396 if (new_rows.size() == 0) {
00397 new_rows.reserve(new_size);
00398 cuts_to_rows(vars, new_cuts, new_rows,
00399 lpres, BCP_Object_FromGenerator, false);
00400 }
00401 for (int i = 0; i < new_size; ++i) {
00402 new_cuts[i]->set_bcpind(-BCP_lp_next_cut_index(*p));
00403 cp.unchecked_push_back(new BCP_lp_waiting_row(new_cuts[i],
00404 new_rows[i]));
00405 }
00406 }
00407 printf("MC: Generated %i cuts from exact solution.\n", new_size);
00408 pool_size = cp.size();
00409
00410 if (par.entry(MC_lp_par::SwitchToSimplex)) {
00411 par.set_entry(MC_lp_par::LpSolver, MC_UseClp);
00412 set_param(BCP_lp_par::MaxCutsAddedPerIteration,
00413 par.entry(MC_lp_par::MaxCutsAddedPerIterSim));
00414 set_param(BCP_lp_par::MaxPresolveIter,
00415 par.entry(MC_lp_par::MaxPresolveIterSim));
00416
00417 getLpProblemPointer()->lp_solver = exact_solver;
00418 getLpProblemPointer()->master_lp = new OsiClpSolverInterface;
00419 getLpProblemPointer()->lp_result->get_results(*exact_solver);
00420 exact_solver = 0;
00421 delete vollp;
00422 vollp = 0;
00423
00424
00425
00426
00427
00428
00429
00430 return BCP_DoNotBranch;
00431 }
00432 }
00433
00434
00435 if (iteration_count > hist_len) {
00436 if (objval < objhist[hist_len-1]) {
00437 objval = objhist[hist_len-1];
00438 }
00439 std::rotate(objhist, objhist+1, objhist+hist_len);
00440 objhist[hist_len-1] = objval;
00441 } else {
00442 if (iteration_count > 1) {
00443 if (objval < objhist[iteration_count-2]) {
00444 objval = objhist[iteration_count-2];
00445 }
00446 }
00447 objhist[iteration_count-1] = objval;
00448 }
00449
00450 if (pool_size != 0 &&
00451 !(tailoff_gap_rel && (tailoff_lb_abs || tailoff_lb_rel))) {
00452 delete exact_solver;
00453 return BCP_DoNotBranch;
00454 }
00455
00456
00457 if (current_level() > par.entry(MC_lp_par::MaxDepth)) {
00458 printf("MC: Maximum depth reached, pruning node.\n");
00459 return BCP_DoNotBranch_Fathomed;
00460 }
00461
00462
00463
00464 if (vollp && !exact_solver) {
00465 exact_solver = solveToOpt(vollp, lpres, vars, cuts, objval);
00466 if (exact_solver == NULL) {
00467 return BCP_DoNotBranch_Fathomed;
00468 }
00469 }
00470
00471
00472
00473 choose_branching_vars(vars, lpres.x(),
00474 par.entry(MC_lp_par::SB_CandidateNum), cands);
00475
00476 if (cands.size() == 0) {
00477
00478
00479
00480
00481 if (local_cut_pool.size() == 0) {
00482 throw BCP_fatal_error("What the heck?! No cuts, integer & tailoff!\n");
00483 }
00484 delete exact_solver;
00485 return BCP_DoNotBranch;
00486 }
00487
00488
00489
00490
00491
00492 if (exact_solver) {
00493 perform_strong_branching(lpres, exact_solver, cands);
00494 delete exact_solver;
00495 }
00496
00497 return BCP_DoBranch;
00498 }
00499
00500
00501
00502 BCP_branching_object_relation
00503 MC_lp::compare_branching_candidates(BCP_presolved_lp_brobj* newobj,
00504 BCP_presolved_lp_brobj* oldobj)
00505 {
00506
00507
00508 return BCP_lp_user::compare_branching_candidates(newobj, oldobj);
00509 }
00510
00511
00512
00513 void
00514 MC_lp::set_actions_for_children(BCP_presolved_lp_brobj* best)
00515 {
00516
00517
00518
00519 if (best_presolved) {
00520 best->swap(*best_presolved);
00521 delete best_presolved;
00522 best_presolved = 0;
00523 }
00524 BCP_lp_user::set_actions_for_children(best);
00525 }
00526
00527
00528
00529
00530 void
00531 MC_lp::perform_strong_branching(const BCP_lp_result& lpres,
00532 OsiSolverInterface* exact_solver,
00533 BCP_vec<BCP_lp_branching_object*>& cands)
00534 {
00535 int i;
00536
00537
00538 exact_solver->markHotStart();
00539
00540
00541 BCP_vec<BCP_lp_branching_object*>::iterator cani;
00542
00543 printf("\n\
00544 MC: Starting strong branching (the objective is transformed!):\
00545 The objective shift is %.4f\n\n", obj_shift);
00546
00547 assert(best_presolved == 0);
00548 int best_fathom_child = 0;
00549 double best_objsum = BCP_DBL_MAX;
00550 BCP_vec<double> tmpobj;
00551
00552 for (cani = cands.begin(); cani != cands.end(); ++cani){
00553
00554 BCP_presolved_lp_brobj* tmp_presolved= new BCP_presolved_lp_brobj(*cani);
00555 const BCP_lp_branching_object* can = *cani;
00556 for (i = 0; i < can->child_num; ++i){
00557 can->apply_child_bd(exact_solver, i);
00558 exact_solver->solveFromHotStart();
00559 tmp_presolved->get_results(*exact_solver, i);
00560 }
00561
00562 tmp_presolved->get_lower_bounds(tmpobj);
00563 for (i = 0; i < can->child_num; ++i) {
00564 tmpobj[i] += obj_shift;
00565 }
00566 tmp_presolved->set_lower_bounds(tmpobj);
00567
00568 const int bvar_ind = (*can->forced_var_pos)[0];
00569 exact_solver->setColBounds(bvar_ind, 0.0, 1.0);
00570 if (get_param(BCP_lp_par::LpVerb_PresolveResult)) {
00571 printf("MC strong branch: Presolving:");
00572 if (get_param(BCP_lp_par::LpVerb_PresolvePositions)) {
00573 can->print_branching_info(exact_solver->getNumCols(), lpres.x(),
00574 exact_solver->getObjCoefficients());
00575 }
00576 for (i = 0; i < can->child_num; ++i) {
00577 const BCP_lp_result& res = tmp_presolved->lpres(i);
00578 const double lb = res.objval();
00579 printf((lb > BCP_DBL_MAX / 10 ? " [%e,%i,%i]" : " [%.4f,%i,%i]"),
00580 lb, res.termcode(), res.iternum());
00581 }
00582 printf("\n");
00583 }
00584
00585
00586
00587 int fathom_child = 0;
00588 double objsum = 0.0;
00589 for (i = 0; i < can->child_num; ++i) {
00590 const BCP_lp_result& res = tmp_presolved->lpres(i);
00591 if (res.termcode() == BCP_ProvenPrimalInf) {
00592 ++fathom_child;
00593 continue;
00594 }
00595 if (res.termcode() == BCP_ProvenOptimal &&
00596 res.objval() + obj_shift + get_param(BCP_lp_par::Granularity) >
00597 upper_bound()) {
00598 ++fathom_child;
00599 continue;
00600 }
00601 objsum += res.objval();
00602 }
00603 if (best_fathom_child < fathom_child) {
00604 std::swap(tmp_presolved, best_presolved);
00605 best_fathom_child = fathom_child;
00606 delete tmp_presolved;
00607 if (best_fathom_child == can->child_num) {
00608 purge_ptr_vector(cands, cani + 1, cands.end());
00609 }
00610 best_objsum = objsum;
00611 continue;
00612 }
00613
00614 if (objsum < best_objsum) {
00615 best_objsum = objsum;
00616 std::swap(tmp_presolved, best_presolved);
00617 }
00618 delete tmp_presolved;
00619 }
00620
00621
00622 exact_solver->unmarkHotStart();
00623
00624
00625 BCP_lp_branching_object* can = best_presolved->candidate();
00626 for (cani = cands.begin(); cani != cands.end(); ++cani) {
00627 if (*cani != can)
00628 delete *cani;
00629 }
00630
00631 cands.clear();
00632 cands.push_back(can);
00633 }