00001
00002
00003 #include <cmath>
00004 #include <algorithm>
00005 #include <numeric>
00006
00007 #include "CoinHelperFunctions.hpp"
00008
00009 #include <OsiVolSolverInterface.hpp>
00010 #include <OsiClpSolverInterface.hpp>
00011
00012 #include "BCP_math.hpp"
00013 #include "BCP_lp.hpp"
00014 #include "BCP_lp_node.hpp"
00015
00016 #include "MC_cut.hpp"
00017 #include "MC_lp.hpp"
00018
00019
00020
00021 static inline void
00022 MC_update_solution(MC_solution*& sol_old, MC_solution*& sol_new)
00023 {
00024 if (!sol_old) {
00025 sol_old = sol_new;
00026 } else {
00027 if (sol_new) {
00028 if (sol_old->objective_value() > sol_new->objective_value()) {
00029 delete sol_old;
00030 sol_old = sol_new;
00031 } else {
00032 delete sol_new;
00033 }
00034 }
00035 }
00036 sol_new = NULL;
00037 }
00038
00039
00040
00041
00042 static inline bool
00043 MC_cycle_row_pair_comp(const std::pair<BCP_cut*, BCP_row*>& p0,
00044 const std::pair<BCP_cut*, BCP_row*>& p1)
00045 {
00046 return MC_cycle_cut_less(p0.first, p1.first);
00047 }
00048
00049
00050
00051 void
00052 MC_lp::unique_cycle_cuts(BCP_vec<BCP_cut*>& new_cuts,
00053 BCP_vec<BCP_row*>& new_rows)
00054 {
00055 size_t i;
00056 const size_t cutnum = new_cuts.size();
00057 if (cutnum != new_rows.size())
00058 throw BCP_fatal_error("Uneven vector sizes in MC_lp_unique_cycle_cuts.\n");
00059
00060 if (cutnum == 0)
00061 return;
00062
00063 std::pair<BCP_cut*, BCP_row*>* cut_row_pairs =
00064 new std::pair<BCP_cut*, BCP_row*>[cutnum];
00065 for (i = 0; i < cutnum; ++i) {
00066 cut_row_pairs[i].first = new_cuts[i];
00067 cut_row_pairs[i].second = new_rows[i];
00068 }
00069
00070 std::sort(cut_row_pairs, cut_row_pairs + cutnum, MC_cycle_row_pair_comp);
00071 int k = 0;
00072 new_cuts[0] = cut_row_pairs[0].first;
00073 new_rows[0] = cut_row_pairs[0].second;
00074 for (i = 1; i < cutnum; ++i) {
00075 if (MC_cycle_cut_equal(new_cuts[k], cut_row_pairs[i].first)) {
00076 delete cut_row_pairs[i].first;
00077 delete cut_row_pairs[i].second;
00078 } else {
00079 new_cuts[++k] = cut_row_pairs[i].first;
00080 new_rows[k] = cut_row_pairs[i].second;
00081 }
00082 }
00083 ++k;
00084
00085 delete[] cut_row_pairs;
00086 new_cuts.erase(new_cuts.entry(k), new_cuts.end());
00087 new_rows.erase(new_rows.entry(k), new_rows.end());
00088 }
00089
00090
00091
00092 void
00093 MC_lp::unpack_module_data(BCP_buffer & buf)
00094 {
00095 par.unpack(buf);
00096 mc.unpack(buf);
00097
00098 mc.create_adj_lists();
00099
00100 hist_len = par.entry(MC_lp_par::TailoffGapRelMinItcount);
00101 hist_len = CoinMax(hist_len, par.entry(MC_lp_par::TailoffLbAbsMinItcount));
00102 hist_len = CoinMax(hist_len, par.entry(MC_lp_par::TailoffLbRelMinItcount));
00103
00104 objhist = hist_len <= 0 ? 0 : new double[hist_len];
00105
00106 if (par.entry(MC_lp_par::LpSolver) & MC_UseVol) {
00107 set_param(BCP_lp_par::MaxCutsAddedPerIteration,
00108 par.entry(MC_lp_par::MaxCutsAddedPerIterVol));
00109 set_param(BCP_lp_par::MaxPresolveIter,
00110 par.entry(MC_lp_par::MaxPresolveIterVol));
00111 } else {
00112 set_param(BCP_lp_par::MaxCutsAddedPerIteration,
00113 par.entry(MC_lp_par::MaxCutsAddedPerIterSim));
00114 set_param(BCP_lp_par::MaxPresolveIter,
00115 par.entry(MC_lp_par::MaxPresolveIterSim));
00116 }
00117 }
00118
00119
00120
00121
00122
00123 OsiSolverInterface *
00124 MC_lp::initialize_solver_interface()
00125 {
00126 if ((par.entry(MC_lp_par::LpSolver) & MC_UseVol) != 0) {
00127 return new OsiVolSolverInterface;
00128 }
00129
00130 if ((par.entry(MC_lp_par::LpSolver) & MC_UseClp) != 0) {
00131 return new OsiClpSolverInterface;
00132 }
00133
00134 throw BCP_fatal_error("MC: No solver is specified!\n");
00135 return 0;
00136 }
00137
00138
00139
00140 void
00141 MC_lp::modify_lp_parameters(OsiSolverInterface* lp, const int changeType,
00142 bool in_strong_branching)
00143 {
00144 if (current_iteration() == 1 &&
00145 ( ((par.entry(MC_lp_par::LpSolver) & MC_UseClp) != 0) )
00146 ) {
00147 started_exact = false;
00148 }
00149 {
00150 OsiVolSolverInterface* vollp = dynamic_cast<OsiVolSolverInterface*>(lp);
00151 if (vollp) {
00152 VOL_parms& vpar = vollp->volprob()->parm;
00153 vpar.lambdainit = par.entry(MC_lp_par::Vol_lambdaInit);
00154 vpar.alphainit = par.entry(MC_lp_par::Vol_alphaInit);
00155 vpar.alphamin = par.entry(MC_lp_par::Vol_alphaMin);
00156 vpar.alphafactor = par.entry(MC_lp_par::Vol_alphaFactor);
00157 vpar.primal_abs_precision = par.entry(MC_lp_par::Vol_primalAbsPrecision);
00158 vpar.gap_abs_precision = par.entry(MC_lp_par::Vol_gapAbsPrecision);
00159 vpar.gap_rel_precision = par.entry(MC_lp_par::Vol_gapRelPrecision);
00160 vpar.granularity = par.entry(MC_lp_par::Vol_granularity);
00161 vpar.minimum_rel_ascent = par.entry(MC_lp_par::Vol_minimumRelAscent);
00162 vpar.ascent_first_check = par.entry(MC_lp_par::Vol_ascentFirstCheck);
00163 vpar.ascent_check_invl = par.entry(MC_lp_par::Vol_ascentCheckInterval);
00164 vpar.maxsgriters = par.entry(MC_lp_par::Vol_maxSubGradientIterations);
00165 vpar.printflag = par.entry(MC_lp_par::Vol_printFlag);
00166 vpar.printinvl = par.entry(MC_lp_par::Vol_printInterval);
00167 vpar.greentestinvl = par.entry(MC_lp_par::Vol_greenTestInterval);
00168 vpar.yellowtestinvl = par.entry(MC_lp_par::Vol_yellowTestInterval);
00169 vpar.redtestinvl = par.entry(MC_lp_par::Vol_redTestInterval);
00170 vpar.alphaint = par.entry(MC_lp_par::Vol_alphaInt);
00171 }
00172 if (in_strong_branching) {
00173
00174 }
00175 }
00176 }
00177
00178
00179
00180 BCP_solution*
00181 MC_lp::test_feasibility(const BCP_lp_result& lp_result,
00182 const BCP_vec<BCP_var*>& vars,
00183 const BCP_vec<BCP_cut*>& cuts)
00184 {
00185 int i, k;
00186
00187 const double etol = par.entry(MC_lp_par::IntegerTolerance);
00188 BCP_vec<double> x(lp_result.x(), vars.size());
00189
00190 for (k = x.size(), i = x.size() - 1; i >= 0; --i) {
00191 const double xx = x[i];
00192 if (xx > 1 - etol) {
00193 --k;
00194 x[i] = 1;
00195 } else if (xx < etol) {
00196 --k;
00197 x[i] = 0;
00198 }
00199 }
00200
00201
00202
00203
00204
00205 if (k == 0) {
00206 BCP_vec<BCP_cut*> new_cuts;
00207 BCP_vec<BCP_row*> new_rows;
00208 const int improve_round = par.entry(MC_lp_par::HeurSwitchImproveRound);
00209 const bool edge_switch = par.entry(MC_lp_par::DoEdgeSwitchHeur);
00210 const int struct_switch = ( par.entry(MC_lp_par::StructureSwitchHeur) &
00211 ((1 << mc.num_structure_type) - 1) );
00212 MC_solution* sol = MC_mst_cutgen(mc, x.begin(), NULL ,
00213 1.0, 0.0,
00214 MC_MstEdgeOrderingPreferExtreme,
00215 improve_round, edge_switch, struct_switch,
00216 .9, 0, new_cuts, new_rows);
00217 purge_ptr_vector(new_rows);
00218 if (new_cuts.size() == 0) {
00219
00220 return sol;
00221 }
00222 purge_ptr_vector(new_cuts);
00223 }
00224
00225 return NULL;
00226 }
00227
00228
00229
00230 BCP_solution*
00231 MC_lp::generate_heuristic_solution(const BCP_lp_result& lpres,
00232 const BCP_vec<BCP_var*>& vars,
00233 const BCP_vec<BCP_cut*>& cuts)
00234 {
00235 return mc_generate_heuristic_solution(lpres.x(), vars, cuts);
00236 }
00237
00238
00239
00240 MC_solution*
00241 MC_lp::mc_generate_heuristic_solution(const double* x,
00242 const BCP_vec<BCP_var*>& vars,
00243 const BCP_vec<BCP_cut*>& cuts)
00244 {
00245 const int m = mc.num_edges;
00246 const int n = mc.num_nodes;
00247
00248 const int heurnum = par.entry(MC_lp_par::MstHeurNum);
00249 const int improve_round = par.entry(MC_lp_par::HeurSwitchImproveRound);
00250 const bool do_edge_switch = par.entry(MC_lp_par::DoEdgeSwitchHeur);
00251 const int struct_switch = ( par.entry(MC_lp_par::StructureSwitchHeur) &
00252 ((1 << mc.num_structure_type) - 1) );
00253 const double maxlambda = par.entry(MC_lp_par::MaxPerturbInMstHeur);
00254
00255 double* w = new double[m];
00256 double * xx = new double[m];
00257
00258 BCP_vec<int> nodesign(n, 1);
00259
00260 MC_solution* sol;
00261 BCP_vec<double> obj;
00262 BCP_vec<char> type;
00263
00264 int i, j;
00265
00266
00267 for (i = 0; i < heurnum; ++i) {
00268 for (j = 0; j < m; ++j)
00269 w[j] = CoinDrand48();
00270 sol = MC_mst_heur(mc, x, w, 1.0, (maxlambda * i) / heurnum,
00271 MC_MstEdgeOrderingPreferExtreme,
00272 improve_round, do_edge_switch, struct_switch);
00273 MC_update_solution(soln, sol);
00274 }
00275 #if 0
00276 for (i = 0; i < heurnum; ++i) {
00277 if (i > 0)
00278 for (j = 0; j < m; ++j)
00279 w[j] = CoinDrand48();
00280 sol = MC_mst_heur(mc, x, w, 1.0, (maxlambda * i) / heurnum,
00281 MC_MstEdgeOrderingPreferZero,
00282 improve_round, do_edge_switch, struct_switch);
00283 MC_update_solution(soln, sol);
00284 }
00285 for (i = 0; i < heurnum; ++i) {
00286 if (i > 0)
00287 for (j = 0; j < m; ++j)
00288 w[j] = CoinDrand48();
00289 sol = MC_mst_heur(mc, x, w, 1.0, (maxlambda * i) / heurnum,
00290 MC_MstEdgeOrderingPreferOne,
00291 improve_round, do_edge_switch, struct_switch);
00292 MC_update_solution(soln, sol);
00293 }
00294 #endif
00295 #if 0
00296
00297 for (i = 0; i < heurnum; ++i) {
00298 for (j = 0; j < m; ++j) {
00299 w[j] = CoinDrand48();
00300 xx[j] = x[j] * (.95 + .1 * CoinDrand48());
00301 }
00302 sol = MC_mst_heur(mc, xx, w, 1.0, (maxlambda * i) / heurnum,
00303 par.entry(MC_lp_par::HeurSwitchImproveRound),
00304 par.entry(MC_lp_par::DoEdgeSwitchHeur),
00305 par.entry(MC_lp_par::DoIsingSquareSwitchHeur));
00306 MC_update_solution(soln, sol);
00307 }
00308 #endif
00309 #if 0
00310
00311
00312 const double * lhs = lpres.lhs();
00313 const int numcuts = cuts.size();
00314 double ratio = 1.0;
00315 for (i = 0; i < numcuts; ++i) {
00316 const double ub = cuts[i]->ub();
00317 if (lhs[i] <= ub)
00318 continue;
00319 if (ub > 0.0) {
00320 ratio = CoinMin(ratio, ub/lhs[i]);
00321 }
00322 }
00323 for (i = 0; i < m; ++i)
00324 xx[i] = x[i]*ratio;
00325 for (i = 0; i < heurnum; ++i) {
00326 for (j = 0; j < m; ++j) {
00327 w[j] = CoinDrand48();
00328 }
00329 sol = MC_mst_heur(mc, xx, w, 1.0, (maxlambda * i) / heurnum,
00330 par.entry(MC_lp_par::HeurSwitchImproveRound),
00331 par.entry(MC_lp_par::DoEdgeSwitchHeur),
00332 par.entry(MC_lp_par::DoIsingSquareSwitchHeur));
00333 MC_update_solution(soln, sol);
00334 }
00335 #endif
00336 delete[] w;
00337 delete[] xx;
00338
00339 sol = soln;
00340 soln = NULL;
00341 return sol;
00342 }
00343
00344
00345
00346 void
00347 MC_lp::pack_feasible_solution(BCP_buffer& buf, const BCP_solution* sol)
00348 {
00349 const MC_solution* mcsol = dynamic_cast<const MC_solution*>(sol);
00350 mcsol->pack(buf);
00351 }
00352
00353
00354
00355 void
00356 MC_lp::cuts_to_rows(const BCP_vec<BCP_var*>& vars,
00357 BCP_vec<BCP_cut*>& cuts,
00358 BCP_vec<BCP_row*>& rows,
00359
00360
00361 const BCP_lp_result& lpres,
00362 BCP_object_origin origin, bool allow_multiple)
00363 {
00364 const int cutnum = cuts.size();
00365 rows.clear();
00366 rows.reserve(cutnum);
00367
00368 const double lb = -BCP_DBL_MAX;
00369
00370 const int varnum = vars.size();
00371 double* coefs = new double[varnum];
00372 int* inds = new int[varnum];
00373 int* iotainds = new int[varnum];
00374 CoinIotaN(iotainds, varnum, 0);
00375
00376 for (int i = 0; i < cutnum; ++i) {
00377 const MC_cycle_cut* cycle_cut =
00378 dynamic_cast<const MC_cycle_cut*>(cuts[i]);
00379 if (cycle_cut) {
00380 const int len = cycle_cut->cycle_len;
00381 const int pos = cycle_cut->pos_edges;
00382 CoinDisjointCopyN(cycle_cut->edges, len, inds);
00383 CoinFillN(coefs, pos, 1.0);
00384 CoinFillN(coefs + pos, len - pos, -1.0);
00385 CoinSort_2(inds, inds + len, coefs);
00386 rows.unchecked_push_back(new BCP_row(inds, inds + len,
00387 coefs, coefs + len,
00388 lb, pos - 1.0));
00389 continue;
00390 }
00391 const MC_explicit_dense_cut* dense_cut =
00392 dynamic_cast<const MC_explicit_dense_cut*>(cuts[i]);
00393 if (dense_cut) {
00394 rows.unchecked_push_back(new BCP_row(iotainds,
00395 iotainds + varnum,
00396 dense_cut->coeffs,
00397 dense_cut->coeffs+varnum,
00398 lb, dense_cut->rhs));
00399 continue;
00400 }
00401 throw BCP_fatal_error("Non-MC_cut to be exanded!\n");
00402 }
00403
00404 if (mc.feas_sol) {
00405 for (int i = 0; i < cutnum; ++i) {
00406 const double lhs = rows[i]->dotProduct(mc.feas_sol->value);
00407 if (cuts[i]->lb() != rows[i]->LowerBound() ||
00408 cuts[i]->ub() != rows[i]->UpperBound()) {
00409 printf("*** What! *** : cut bounds are not the same: %i %f %f %f %f\n",
00410 i, cuts[i]->lb(), rows[i]->LowerBound(),
00411 cuts[i]->ub(), rows[i]->UpperBound());
00412 }
00413 if (lhs < cuts[i]->lb() || lhs > cuts[i]->ub()) {
00414 printf("*** Violation *** : index: %i lhs: %f lb: %f ub: %f\n",
00415 i, lhs, cuts[i]->lb(), cuts[i]->ub());
00416 }
00417 }
00418 }
00419
00420 delete[] iotainds;
00421 delete[] inds;
00422 delete[] coefs;
00423 }
00424
00425
00426
00427 BCP_object_compare_result
00428 MC_lp::compare_cuts(const BCP_cut* c0, const BCP_cut* c1)
00429 {
00430 return MC_cycle_cut_equal(c0, c1) ? BCP_ObjsAreSame : BCP_DifferentObjs;
00431 }
00432
00433
00434
00435
00436 void
00437 MC_lp::generate_mst_cuts(const double* x, const double* lhs,
00438 const double objval,
00439 const BCP_vec<BCP_var*>& vars,
00440 const BCP_vec<BCP_cut*>& cuts,
00441 BCP_vec<BCP_cut*>& new_cuts,
00442 BCP_vec<BCP_row*>& new_rows)
00443 {
00444 const double max_lp_viol = getMaxLpViol();
00445 const int heurnum = par.entry(MC_lp_par::CycleCutHeurNum);
00446 const double minviol =
00447 CoinMax(max_lp_viol, par.entry(MC_lp_par::MinMstCycleCutViolation));
00448 const double maxlambda = par.entry(MC_lp_par::MaxPerturbInMstCycleCutGen);
00449
00450 const int heurswitchround = par.entry(MC_lp_par::HeurSwitchImproveRound);
00451 const bool edge_switch = par.entry(MC_lp_par::DoEdgeSwitchHeur);
00452 const int struct_switch = ( par.entry(MC_lp_par::StructureSwitchHeur) &
00453 ((1 << mc.num_structure_type) - 1) );
00454
00455 const int numvars = vars.size();
00456 double* w = new double[numvars];
00457 double* xx = new double[numvars];
00458 const int max_cycle_num = par.entry(MC_lp_par::MaxCycleCutNum);
00459
00460
00461 int i, j;
00462 MC_solution* sol = NULL;
00463 for (i = 0; i < heurnum; ++i) {
00464 if (i > 0)
00465 for (j = 0; j < numvars; ++j)
00466 w[j] = CoinDrand48();
00467 sol = MC_mst_cutgen(mc, x, w, 1.0, (maxlambda*i)/heurnum,
00468 MC_MstEdgeOrderingPreferExtreme,
00469 heurswitchround, edge_switch, struct_switch,
00470 minviol, new_cuts.size() + max_cycle_num,
00471 new_cuts, new_rows);
00472 MC_update_solution(soln, sol);
00473 }
00474 #if 0
00475 for (i = 0; i < heurnum; ++i) {
00476 if (i > 0)
00477 for (j = 0; j < numvars; ++j)
00478 w[j] = CoinDrand48();
00479 sol = MC_mst_cutgen(mc, x, w, 1.0, (maxlambda*i)/heurnum,
00480 MC_MstEdgeOrderingPreferOne, minviol,
00481 heurswitchround, edge_switch, struct_switch,
00482 minviol, new_cuts.size() + max_cycle_num,
00483 new_cuts, new_rows);
00484 delete sol;
00485 }
00486 for (i = 0; i < heurnum; ++i) {
00487 if (i > 0)
00488 for (j = 0; j < numvars; ++j)
00489 w[j] = CoinDrand48();
00490 sol = MC_mst_cutgen(mc, x, w, 1.0, (maxlambda*i)/heurnum,
00491 MC_MstEdgeOrderingPreferExtreme, minviol,
00492 heurswitchround, edge_switch, struct_switch,
00493 minviol, new_cuts.size() + max_cycle_num,
00494 new_cuts, new_rows);
00495 delete sol;
00496 }
00497 const int cutnum_add = new_cuts.size() - cutnum;
00498 printf("MC: cycle cuts (add): %i\n", cutnum_add);
00499 #endif
00500 #if 0
00501
00502 for (i = 0; i < numvars; ++i)
00503 xx[i] = x[i]*(.95 + .1 * CoinDrand48());
00504 for (i = 0; i < heurnum; ++i) {
00505 if (i > 0) {
00506 for (j = 0; j < numvars; ++j)
00507 w[j] = CoinDrand48();
00508 }
00509 sol = MC_mst_cutgen(mc, xx, w, 1.0, (maxlambda*i)/heurnum, minviol,
00510 heurswitchround, new_cuts.size() + max_cycle_num,
00511 new_cuts, new_rows);
00512 printf("MC: sol value: %.0f\n", sol->objective_value());
00513
00514 MC_update_solution(soln, sol);
00515 }
00516 const int cutnum_mult_add = new_cuts.size() - cutnum - cutnum_add;
00517 printf("MC: cycle cuts (mult & add): %i\n", cutnum_mult_add);
00518 #endif
00519
00520
00521
00522 {
00523 #if 0
00524 const double * lhs = lpres.lhs();
00525 const int numcuts = cuts.size();
00526 double ratio = 1.0;
00527 int cnt = 0;
00528 double maxviol = 0.0;
00529 for (i = 0; i < numcuts; ++i) {
00530 const double ub = cuts[i]->ub();
00531 if (lhs[i] <= ub)
00532 continue;
00533 if (ub == 0.0) {
00534 ++cnt;
00535 maxviol = CoinMax(maxviol, lhs[i] - ub);
00536 } else {
00537 ratio = CoinMin(ratio, ub/lhs[i]);
00538 }
00539 }
00540 printf("MC: primal shrinking: %.6f , %i cuts left violated (%.6f)\n",
00541 ratio, cnt, maxviol);
00542 for (i = 0; i < numvars; ++i)
00543 xx[i] = x[i]*ratio;
00544 for (i = 0; i < heurnum; ++i) {
00545 if (i > 0) {
00546 for (j = 0; j < numvars; ++j)
00547 w[j] = CoinDrand48();
00548 }
00549 sol = MC_mst_cutgen(mc, xx, w, 1.0, (maxlambda*i)/heurnum, minviol,
00550 heurswitchround, new_cuts.size() + max_cycle_num,
00551 new_cuts, new_rows);
00552 printf("MC: sol value: %.0f\n", sol->objective_value());
00553 MC_update_solution(soln, sol);
00554 }
00555 const int cutnum_shrink_add =
00556 new_cuts.size() - cutnum - cutnum_add - cutnum_mult_add;
00557 printf("MC: cycle cuts (shrink & add): %i\n", cutnum_shrink_add);
00558 #endif
00559 }
00560 delete[] w;
00561 delete[] xx;
00562 }
00563
00564 void
00565 MC_lp::generate_sp_cuts(const double* x, const double* lhs,
00566 const double objval,
00567 const BCP_vec<BCP_var*>& vars,
00568 const BCP_vec<BCP_cut*>& cuts,
00569 BCP_vec<BCP_cut*>& new_cuts,
00570 BCP_vec<BCP_row*>& new_rows)
00571 {
00572 const double max_lp_viol = getMaxLpViol();
00573 const bool all_sp_cuts = par.entry(MC_lp_par::ReportAllSPCycleCuts);
00574 const double minviol =
00575 CoinMax(max_lp_viol, par.entry(MC_lp_par::MinSPCycleCutViolation));
00576 MC_generate_shortest_path_cycles(mc, x, all_sp_cuts, minviol,
00577 new_cuts, new_rows);
00578 }
00579
00580 void
00581 MC_lp::generate_cuts_in_lp(const BCP_lp_result& lpres,
00582 const BCP_vec<BCP_var*>& vars,
00583 const BCP_vec<BCP_cut*>& cuts,
00584 BCP_vec<BCP_cut*>& new_cuts,
00585 BCP_vec<BCP_row*>& new_rows)
00586 {
00587 #if 0
00588 const double* x = lpres.x();
00589 double* xx = new double[vars.size()];
00590 for (int i = vars.size() - 1; i >= 0; --i) {
00591 if (x[i] < 0.0)
00592 xx[i] = 0.0;
00593 else if (x[i] > 1.0)
00594 xx[i] = 1.0;
00595 else
00596 xx[i] = x[i];
00597 }
00598 generate_cuts_in_lp(xx, lpres.lhs(), lpres.objval(),
00599 vars, cuts, new_cuts, new_rows);
00600 delete[] xx;
00601 #else
00602 generate_cuts_in_lp(lpres.x(), lpres.lhs(), lpres.objval(),
00603 vars, cuts, new_cuts, new_rows);
00604 #endif
00605 }
00606
00607 double
00608 MC_lp::getMaxLpViol()
00609 {
00610 double max_lp_viol = 0.0;
00611
00612
00613
00614 OsiVolSolverInterface* vollp =
00615 dynamic_cast<OsiVolSolverInterface*>(getLpProblemPointer()->lp_solver);
00616 if (vollp) {
00617 const double * lhs = getLpProblemPointer()->lp_result->lhs();
00618 const double * rhs = vollp->getRightHandSide();
00619 for (int k = vollp->getNumRows() - 1; k >= 0; --k)
00620 max_lp_viol = CoinMax(max_lp_viol, lhs[k] - rhs[k]);
00621 }
00622 max_lp_viol += 0.001;
00623 return max_lp_viol;
00624 }
00625
00626 void
00627 MC_lp::generate_cuts_in_lp(const double* x, const double* lhs,
00628 const double objval,
00629 const BCP_vec<BCP_var*>& vars,
00630 const BCP_vec<BCP_cut*>& cuts,
00631 BCP_vec<BCP_cut*>& new_cuts,
00632 BCP_vec<BCP_row*>& new_rows)
00633 {
00634 bool tailoff_gap_rel, tailoff_lb_abs, tailoff_lb_rel;
00635 tailoff_test(tailoff_gap_rel, tailoff_lb_abs, tailoff_lb_rel, objval);
00636
00637 const double max_lp_viol = getMaxLpViol();
00638
00639 int i, cutnum = new_cuts.size();
00640
00641 if (mc.ising_four_cycles || mc.ising_triangles) {
00642 const int grid = static_cast<int>(sqrt(mc.num_nodes + 1.0));
00643 const int grid_nodes = grid*grid;
00644 const double minviol =
00645 CoinMax(max_lp_viol, par.entry(MC_lp_par::MinIsingCutViolation));
00646
00647 if (mc.ising_four_cycles) {
00648 MC_test_ising_four_cycles(grid_nodes, mc.ising_four_cycles,
00649 x, minviol, new_cuts, new_rows);
00650 const int ising_four_cycle_num = new_cuts.size() - cutnum;
00651 printf("MC: ising four cycles: %i \n", ising_four_cycle_num);
00652 cutnum = new_cuts.size();
00653 }
00654
00655 if (mc.ising_triangles) {
00656 MC_test_ising_triangles(2*grid_nodes, mc.ising_triangles,
00657 x, 0.9, new_cuts, new_rows);
00658 const int ising_triangle_num = new_cuts.size() - cutnum;
00659 printf("MC: ising triangles: %i \n", ising_triangle_num);
00660 cutnum = new_cuts.size();
00661 }
00662 }
00663
00664 const int gen_mst_cycle = par.entry(MC_lp_par::MstCycleCutGeneration);
00665 #if 0
00666 const bool mst = gen_mst_cycle == MC_AlwaysGenerateMstCycleCuts;
00667 #else
00668 const bool mst =
00669 gen_mst_cycle == MC_AlwaysGenerateMstCycleCuts ||
00670 (gen_mst_cycle == MC_GenerateMstCycleCutsAsLastResort &&
00671 (cutnum < 50 || tailoff_gap_rel || tailoff_lb_abs || tailoff_lb_rel));
00672 #endif
00673
00674 if (mst) {
00675 generate_mst_cuts(x, lhs, objval, vars, cuts, new_cuts, new_rows);
00676 const int cycle_num = new_cuts.size() - cutnum;
00677 unique_cycle_cuts(new_cuts, new_rows);
00678 printf("MC: cycle cuts: %i (%i before removing duplicates)\n",
00679 static_cast<int>(new_cuts.size()) - cutnum, cycle_num);
00680 cutnum = new_cuts.size();
00681 }
00682
00683 const MC_SPCycleCutGen gen_sp_cycle =
00684 static_cast<MC_SPCycleCutGen>(par.entry(MC_lp_par::SPCycleCutGeneration));
00685 #if 0
00686 const bool sp_cuts = gen_sp_cycle == MC_AlwaysGenerateSPCycleCuts;
00687 #else
00688 const bool sp_cuts =
00689 gen_sp_cycle == MC_AlwaysGenerateSPCycleCuts ||
00690 (gen_sp_cycle == MC_GenerateSPCycleCutsAsLastResort &&
00691 (cutnum < 50 || tailoff_gap_rel || tailoff_lb_abs || tailoff_lb_rel));
00692 #endif
00693
00694 if (sp_cuts) {
00695 generate_sp_cuts(x, lhs, objval, vars, cuts, new_cuts, new_rows);
00696 const int sp_cycle_num = new_cuts.size() - cutnum;
00697 unique_cycle_cuts(new_cuts, new_rows);
00698 printf("MC: SP based cycle cuts: %i (%i before removing duplicates)\n",
00699 static_cast<int>(new_cuts.size()) - cutnum, sp_cycle_num);
00700 cutnum = new_cuts.size();
00701 }
00702
00703 if (mc.feas_sol) {
00704 for (i = 0; i < cutnum; ++i) {
00705 const double lhs = new_rows[i]->dotProduct(mc.feas_sol->value);
00706 if (new_cuts[i]->lb() != new_rows[i]->LowerBound() ||
00707 new_cuts[i]->ub() != new_rows[i]->UpperBound()) {
00708 printf("*** What! *** : cut bounds are not the same: %i %f %f %f %f\n",
00709 i, new_cuts[i]->lb(), new_rows[i]->LowerBound(),
00710 new_cuts[i]->ub(), new_rows[i]->UpperBound());
00711 }
00712 if (lhs < new_cuts[i]->lb() || lhs > new_cuts[i]->ub()) {
00713 printf("*** Violation *** : index: %i lhs: %f lb: %f ub: %f\n",
00714 i, lhs, new_cuts[i]->lb(), new_cuts[i]->ub());
00715 }
00716 }
00717 }
00718 }
00719
00720
00721
00722 void
00723 MC_lp::logical_fixing(const BCP_lp_result& lpres,
00724 const BCP_vec<BCP_var*>& vars,
00725 const BCP_vec<BCP_cut*>& cuts,
00726 const BCP_vec<BCP_obj_status>& var_status,
00727 const BCP_vec<BCP_obj_status>& cut_status,
00728 const int var_bound_changes_since_logical_fixing,
00729 BCP_vec<int>& changed_pos, BCP_vec<double>& new_bd)
00730 {
00731
00732
00733
00734 if (var_bound_changes_since_logical_fixing == 0)
00735 return;
00736 #if 1
00737 const int n = mc.num_nodes;
00738 const int m = mc.num_edges;
00739 const MC_graph_edge* edges = mc.edges;
00740
00741 int * first_on_chain = new int[n];
00742 int * last_on_chain = new int[n];
00743 int * next_on_chain = new int[n];
00744 int * size_of_chain = new int[n];
00745 int * sig = new int[n];
00746
00747 CoinIotaN(first_on_chain, n, 0);
00748 CoinIotaN(last_on_chain, n, 0);
00749 CoinFillN(next_on_chain, n, -1);
00750 CoinFillN(size_of_chain, n, 1);
00751 CoinFillN(sig, n, 1);
00752
00753 int tree_size = 0;
00754
00755 int label_s = -1;
00756 int label_l = -1;
00757 for (int k = 0; k < m; ++k) {
00758 if (vars[k]->lb() != vars[k]->ub())
00759 continue;
00760 const int i = edges[k].tail;
00761 const int j = edges[k].head;
00762 const int label_i = first_on_chain[i];
00763 const int label_j = first_on_chain[j];
00764 if (label_i == label_j)
00765 continue;
00766 if (size_of_chain[label_i] > size_of_chain[label_j]) {
00767 label_s = label_j;
00768 label_l = label_i;
00769 } else {
00770 label_s = label_i;
00771 label_l = label_j;
00772 }
00773 ++tree_size;
00774 for (int l = label_s; l != -1; l = next_on_chain[l]) {
00775 first_on_chain[l] = label_l;
00776 }
00777 if ((vars[k]->lb() == 0.0 && sig[i] != sig[j]) ||
00778 (vars[k]->lb() == 1.0 && sig[i] == sig[j])) {
00779
00780 for (int l = label_s; l != -1; l = next_on_chain[l]) {
00781 sig[l] *= -1;
00782 }
00783 }
00784 next_on_chain[last_on_chain[label_l]] = label_s;
00785 last_on_chain[label_l] = last_on_chain[label_s];
00786 size_of_chain[label_l] += size_of_chain[label_s];
00787 }
00788
00789 delete[] last_on_chain;
00790 delete[] next_on_chain;
00791 delete[] size_of_chain;
00792
00793
00794
00795
00796
00797 const int * component = first_on_chain;
00798 for (int k = 0; k < m; ++k) {
00799 if (vars[k]->lb() == vars[k]->ub())
00800 continue;
00801 const int i = edges[k].tail;
00802 const int j = edges[k].head;
00803 if (component[i] == component[j]) {
00804
00805 changed_pos.push_back(k);
00806 if (sig[i] == sig[j]) {
00807 new_bd.push_back(0.0);
00808 new_bd.push_back(0.0);
00809 } else {
00810 new_bd.push_back(1.0);
00811 new_bd.push_back(1.0);
00812 }
00813 }
00814 }
00815 printf("MC: logical fixing: # of components: %i, fixed %i variables.\n",
00816 n - tree_size, static_cast<int>(changed_pos.size()));
00817 delete[] first_on_chain;
00818 delete[] sig;
00819 #endif
00820 }