00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "BonCbc.hpp"
00012 #include "BonBabInfos.hpp"
00013 #include "CglCutGenerator.hpp"
00014
00015 #include "CouenneCutGenerator.hpp"
00016 #include "CouenneProblem.hpp"
00017 #include "CouenneProblemElem.hpp"
00018 #include "CouenneExprVar.hpp"
00019 #include "CouenneInfeasCut.hpp"
00020
00021 #include "CouenneRecordBestSol.hpp"
00022
00023 #ifdef COIN_HAS_NTY
00024 #include "Nauty.h"
00025 #endif
00026
00027 using namespace Ipopt;
00028
00029 namespace Couenne {
00030
00031 #define Couenne_large_bound2 9.99e12
00032
00033
00034 bool isOptimumCut (const CouNumber *opt, OsiCuts &cs, CouenneProblem *p);
00035
00036
00037
00038 void fictitiousBound (OsiCuts &cs,
00039 CouenneProblem *p,
00040 bool action) {
00041
00042
00043 const CouNumber large_tol = (Couenne_large_bound2 / 1e6);
00044
00045
00046
00047 int ind_obj = p -> Obj (0) -> Body () -> Index ();
00048
00049 if (ind_obj < 0) return;
00050
00051
00052
00053
00054
00055 if (action)
00056
00057 {if (p -> Lb (ind_obj) < - Couenne_large_bound2) p -> Lb (ind_obj) = - Couenne_large_bound2;}
00058
00059 else
00060
00061
00062 {if (fabs (p->Lb(ind_obj)+Couenne_large_bound2)<large_tol) p->Lb(ind_obj) =-COUENNE_INFINITY;}
00063 }
00064
00065
00066
00067 void sparse2dense (int ncols, t_chg_bounds *chg_bds, int *&changed, int &nchanged) {
00068
00069
00070
00071 changed = (int *) realloc (changed, ncols * sizeof (int));
00072 nchanged = 0;
00073
00074 for (register int i=ncols, j=0; i--; j++, chg_bds++)
00075 if (chg_bds -> lower() != t_chg_bounds::UNCHANGED ||
00076 chg_bds -> upper() != t_chg_bounds::UNCHANGED ) {
00077 *changed++ = j;
00078 nchanged++;
00079 }
00080
00081 changed -= nchanged;
00082
00083 }
00084
00085
00087 void updateBranchInfo (const OsiSolverInterface &si, CouenneProblem *p,
00088 t_chg_bounds *chg, const CglTreeInfo &info);
00089
00091
00092 void CouenneCutGenerator::generateCuts (const OsiSolverInterface &si,
00093 OsiCuts &cs,
00094 const CglTreeInfo info)
00095 #if CGL_VERSION_MAJOR == 0 && CGL_VERSION_MINOR <= 57
00096 const
00097 #endif
00098 {
00099
00100
00101
00102
00103
00104 if (isWiped (cs) ||
00105 (CoinCpuTime () > problem_ -> getMaxCpuTime ()))
00106 return;
00107
00108 #ifdef FM_TRACE_OPTSOL
00109 double currCutOff = problem_->getCutOff();
00110 double bestVal = 1e50;
00111 CouenneRecordBestSol *rs = problem_->getRecordBestSol();
00112 if(rs->getHasSol()) {
00113 bestVal = rs->getVal();
00114 }
00115 if(currCutOff > bestVal) {
00116
00117 problem_ -> setCutOff (bestVal);
00118
00119 int indObj = problem_->Obj(0)->Body()->Index();
00120
00121 if (indObj >= 0) {
00122 OsiColCut *objCut = new OsiColCut;
00123 objCut->setUbs(1, &indObj, &bestVal);
00124 cs.insert(objCut);
00125 delete objCut;
00126 }
00127 }
00128 #endif
00129
00130 #ifdef FM_PRINT_INFO
00131 if((BabPtr_ != NULL) && (info.level >= 0) && (info.pass == 0) &&
00132 (BabPtr_->model().getNodeCount() > lastPrintLine)) {
00133 printLineInfo();
00134 lastPrintLine += 1;
00135 }
00136 #endif
00137
00138 const int infeasible = 1;
00139
00140 int nInitCuts = cs.sizeRowCuts ();
00141
00142 CouNumber
00143 *&realOpt = problem_ -> bestSol (),
00144 *saveOptimum = realOpt;
00145
00146 if (!firstcall_ && realOpt) {
00147
00148
00149
00150
00151 CouNumber *opt = realOpt;
00152
00153 const CouNumber
00154 *sol = si.getColSolution (),
00155 *lb = si.getColLower (),
00156 *ub = si.getColUpper ();
00157
00158 int objind = problem_ -> Obj (0) -> Body () -> Index ();
00159
00160 for (int j=0, i=problem_ -> nVars (); i--; j++, opt++, lb++, ub++)
00161 if ((j != objind) &&
00162 ((*opt < *lb - COUENNE_EPS * (1 + CoinMin (fabs (*opt), fabs (*lb)))) ||
00163 (*opt > *ub + COUENNE_EPS * (1 + CoinMin (fabs (*opt), fabs (*ub)))))) {
00164
00165 jnlst_ -> Printf (J_VECTOR, J_CONVEXIFYING,
00166 "out of bounds, ignore x%d = %g [%g,%g] opt = %g\n",
00167 problem_ -> nVars () - i - 1, *sol, *lb, *ub, *opt);
00168
00169
00170
00171 realOpt = NULL;
00172 break;
00173 }
00174 }
00175
00176
00177
00178
00179
00180
00181
00182 jnlst_ -> Printf (J_DETAILED, J_CONVEXIFYING,
00183 "generateCuts: level = %d, pass = %d, intree = %d\n",
00184 info.level, info.pass, info.inTree);
00185
00186 Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (si.getAuxiliaryInfo ());
00187
00188 if (babInfo)
00189 babInfo -> setFeasibleNode ();
00190
00191 double now = CoinCpuTime ();
00192 int ncols = problem_ -> nVars ();
00193
00194
00195
00196
00197
00198 t_chg_bounds *chg_bds = new t_chg_bounds [ncols];
00199
00200
00201
00202
00203
00204
00205
00206 problem_ -> installCutOff ();
00207
00208 if (firstcall_) {
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218 for (int i=0; i < ncols; i++)
00219 if (problem_ -> Var (i) -> Multiplicity () > 0) {
00220 chg_bds [i].setLower (t_chg_bounds::CHANGED);
00221 chg_bds [i].setUpper (t_chg_bounds::CHANGED);
00222 }
00223
00224
00225
00226 if (problem_ -> doFBBT () &&
00227 (! (problem_ -> boundTightening (chg_bds, babInfo))))
00228 jnlst_ -> Printf (J_STRONGWARNING, J_CONVEXIFYING,
00229 "Couenne: WARNING, first convexification is infeasible\n");
00230
00231
00232
00233
00234
00235 int nnlc = problem_ -> nCons ();
00236
00237 for (int i=0; i<nnlc; i++) {
00238
00239 if (CoinCpuTime () > problem_ -> getMaxCpuTime ())
00240 break;
00241
00242
00243 CouenneConstraint *con = problem_ -> Con (i);
00244
00245
00246 int objindex = con -> Body () -> Index ();
00247
00248 if ((objindex >= 0) &&
00249 ((con -> Body () -> Type () == AUX) ||
00250 (con -> Body () -> Type () == VAR))) {
00251
00252
00253 exprVar *conaux = problem_ -> Var (objindex);
00254
00255 if (conaux &&
00256 (conaux -> Type () == AUX) &&
00257 (conaux -> Image ()) &&
00258 (conaux -> Image () -> Linearity () <= LINEAR)) {
00259
00260
00261
00262
00263 double
00264 lb = (*(con -> Lb ())) (),
00265 ub = (*(con -> Ub ())) ();
00266
00267 OsiColCut newBound;
00268 if (lb > -COUENNE_INFINITY) newBound.setLbs (1, &objindex, &lb);
00269 if (ub < COUENNE_INFINITY) newBound.setUbs (1, &objindex, &ub);
00270
00271 cs.insert (newBound);
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290 }
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 } else {
00305
00306
00307 assert (false);
00308 }
00309 }
00310
00311 if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_CONVEXIFYING)) {
00312 if (cs.sizeRowCuts ()) {
00313 jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne: %d constraint row cuts\n",
00314 cs.sizeRowCuts ());
00315 for (int i=0; i<cs.sizeRowCuts (); i++)
00316 cs.rowCutPtr (i) -> print ();
00317 }
00318 if (cs.sizeColCuts ()) {
00319 jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne: %d constraint col cuts\n",
00320 cs.sizeColCuts ());
00321 for (int i=0; i<cs.sizeColCuts (); i++)
00322 cs.colCutPtr (i) -> print ();
00323 }
00324 }
00325 } else {
00326
00327
00328 int indobj = problem_ -> Obj (0) -> Body () -> Index ();
00329
00330
00331 problem_ -> domain () -> push (&si, &cs);
00332
00333 if (indobj >= 0) {
00334
00335
00336
00337 double lp_bound = problem_ -> domain () -> x (indobj);
00338
00339
00340 {if (lp_bound > problem_ -> Lb (indobj)) problem_ -> Lb (indobj) = lp_bound;}
00341
00342 }
00343
00344 updateBranchInfo (si, problem_, chg_bds, info);
00345 }
00346
00347
00348 for (int i = problem_ -> nCons (); i--;) {
00349
00350
00351 CouenneConstraint *con = problem_ -> Con (i);
00352
00353
00354 int objindex = con -> Body () -> Index ();
00355
00356 if ((objindex >= 0) &&
00357 ((con -> Body () -> Type () == AUX) ||
00358 (con -> Body () -> Type () == VAR))) {
00359
00360
00361 CouNumber
00362 l = con -> Lb () -> Value (),
00363 u = con -> Ub () -> Value ();
00364
00365
00366 problem_ -> Lb (objindex) = CoinMax (l, problem_ -> Lb (objindex));
00367 problem_ -> Ub (objindex) = CoinMin (u, problem_ -> Ub (objindex));
00368 }
00369 }
00370
00371 problem_ -> installCutOff ();
00372
00373 fictitiousBound (cs, problem_, false);
00374
00375 int *changed = NULL, nchanged;
00376
00377
00378
00379
00380
00381
00382
00383 try {
00384
00385
00386
00387
00388
00389 #ifdef COIN_HAS_NTY
00390
00391
00392
00393 if (problem_ -> orbitalBranching ())
00394 problem_ -> Compute_Symmetry ();
00395 #endif
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 if (!firstcall_ &&
00408 problem_ -> doRCBT () &&
00409 problem_ -> redCostBT (&si, chg_bds) &&
00410 !(problem_ -> btCore (chg_bds)))
00411 throw infeasible;
00412
00413
00414 if (problem_ -> doFBBT () &&
00415
00416 (! (problem_ -> boundTightening (chg_bds, babInfo))))
00417 throw infeasible;
00418
00419
00420 if (!firstcall_ &&
00421 problem_ -> obbt (this, si, cs, info, babInfo, chg_bds) < 0)
00422 throw infeasible;
00423
00424
00425
00426 if ((problem_ -> doFBBT () ||
00427 problem_ -> doOBBT () ||
00428 problem_ -> doABT ()) &&
00429 (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING))) {
00430
00431 jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== after bt =============\n");
00432 for (int i = 0; i < problem_ -> nVars (); i++)
00433 if (problem_ -> Var (i) -> Multiplicity () > 0)
00434 jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
00435 problem_ -> X (i), problem_ -> Lb (i), problem_ -> Ub (i));
00436 jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
00437 }
00438
00439
00440
00441 #ifdef COIN_HAS_NTY
00442
00443
00444
00445
00446 if (problem_ -> orbitalBranching () && !firstcall_) {
00447
00448 CouNumber
00449 *lb = problem_ -> Lb (),
00450 *ub = problem_ -> Ub ();
00451
00452 std::vector<std::vector<int> > *new_orbits = problem_ -> getNtyInfo () -> getOrbits();
00453
00454 for (int i=0, ii = problem_ -> getNtyInfo () -> getNumOrbits (); ii--; i++){
00455
00456 CouNumber
00457 ll = -COUENNE_INFINITY,
00458 uu = COUENNE_INFINITY;
00459
00460 std::vector <int> orbit = (*new_orbits)[i];
00461
00462 if (orbit.size () <= 1)
00463 continue;
00464
00465 if (jnlst_ -> ProduceOutput (J_VECTOR, J_BOUNDTIGHTENING)) {
00466 printf ("orbit bounds: "); fflush (stdout);
00467 for(int j = 0; j < orbit.size (); j++) {
00468 printf ("x_%d [%g,%g] ", orbit[j], lb [orbit [j]], ub [orbit [j]]);
00469 fflush (stdout);
00470 }
00471 printf ("\n");
00472 }
00473
00474 for (int j = 0; j < orbit.size (); j++) {
00475
00476 int indOrb = orbit [j];
00477
00478 if (indOrb < problem_ -> nVars ()) {
00479
00480 if (lb [indOrb] > ll) ll = lb [indOrb];
00481 if (ub [indOrb] < uu) uu = ub [indOrb];
00482 }
00483 }
00484
00485 jnlst_ -> Printf (J_VECTOR, J_BOUNDTIGHTENING,
00486 " --> new common lower bounds: [%g,--]\n", ll);
00487
00488 for(int j = 0; j < orbit.size (); j++) {
00489
00490 int indOrb = orbit [j];
00491
00492 if (indOrb < problem_ -> nVars ()){
00493
00494 lb [indOrb] = ll;
00495 ub [indOrb] = uu;
00496 }
00497 }
00498 }
00499
00500 delete new_orbits;
00501 }
00502
00503 #endif
00504
00505
00506
00507 sparse2dense (ncols, chg_bds, changed, nchanged);
00508
00509 double *nlpSol = NULL;
00510
00511
00512
00513 if (true) {
00514
00515 if (babInfo)
00516 nlpSol = const_cast <double *> (babInfo -> nlpSolution ());
00517
00518
00519
00520 int logAbtLev = problem_ -> logAbtLev ();
00521
00522 if (problem_ -> doABT () &&
00523 ((logAbtLev != 0) ||
00524 (info.level == 0)) &&
00525 (info.pass == 0) &&
00526 ((logAbtLev < 0) ||
00527 (info.level <= logAbtLev) ||
00528 (CoinDrand48 () <
00529 pow (2., (double) logAbtLev - (info.level + 1))))) {
00530
00531 jnlst_ -> Printf(J_VECTOR, J_BOUNDTIGHTENING," performing ABT\n");
00532 if (! (problem_ -> aggressiveBT (nlp_, chg_bds, info, babInfo)))
00533 throw infeasible;
00534
00535 sparse2dense (ncols, chg_bds, changed, nchanged);
00536 }
00537
00538
00539
00540
00541
00542
00543
00544
00545 bool save_av = addviolated_;
00546 addviolated_ = false;
00547
00548
00549 problem_ -> domain () -> push
00550 (problem_ -> nVars (),
00551 problem_ -> domain () -> x (),
00552 problem_ -> domain () -> lb (),
00553 problem_ -> domain () -> ub (), false);
00554
00555
00556 if (nlpSol) {
00557 CoinCopyN (nlpSol, problem_ -> nOrigVars (), problem_ -> domain () -> x ());
00558
00559
00560 problem_ -> getAuxs (problem_ -> domain () -> x ());
00561 }
00562
00563 if (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING)) {
00564 jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== genrowcuts on NLP =============\n");
00565 for (int i = 0; i < problem_ -> nVars (); i++)
00566 if (problem_ -> Var (i) -> Multiplicity () > 0)
00567 jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
00568 problem_ -> X (i),
00569 problem_ -> Lb (i),
00570 problem_ -> Ub (i));
00571 jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
00572 }
00573
00574 problem_ -> domain () -> current () -> isNlp () = true;
00575 genRowCuts (si, cs, nchanged, changed, chg_bds);
00576
00577 problem_ -> domain () -> pop ();
00578
00579 addviolated_ = save_av;
00580
00581
00582 if (babInfo)
00583 babInfo -> setHasNlpSolution (false);
00584
00585 } else {
00586
00587 if (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING)) {
00588 jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== genrowcuts on LP =============\n");
00589 for (int i = 0; i < problem_ -> nVars (); i++)
00590 if (problem_ -> Var (i) -> Multiplicity () > 0)
00591 jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
00592 problem_ -> X (i),
00593 problem_ -> Lb (i),
00594 problem_ -> Ub (i));
00595 jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
00596 }
00597
00598 genRowCuts (si, cs, nchanged, changed, chg_bds);
00599 }
00600
00601
00602 if (nchanged)
00603 genColCuts (si, cs, nchanged, changed);
00604
00605 if (firstcall_ && (cs.sizeRowCuts () >= 1))
00606 jnlst_->Printf(J_ITERSUMMARY, J_CONVEXIFYING,
00607 "Couenne: %d initial row cuts\n", cs.sizeRowCuts ());
00608
00609 if (realOpt &&
00610 isOptimumCut (realOpt, cs, problem_))
00611 jnlst_->Printf(J_ITERSUMMARY, J_CONVEXIFYING,
00612 "Warning: Optimal solution was cut\n");
00613 }
00614
00615 catch (int exception) {
00616
00617 if ((exception == infeasible) && (!firstcall_)) {
00618
00619 jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,
00620 "Couenne: Infeasible node\n");
00621
00622 WipeMakeInfeas (cs);
00623 }
00624
00625 if (babInfo)
00626 babInfo -> setInfeasibleNode ();
00627 }
00628
00629 delete [] chg_bds;
00630
00631 if (changed)
00632 free (changed);
00633
00634 if (firstcall_) {
00635
00636 jnlst_ -> Printf (J_SUMMARY, J_CONVEXIFYING,
00637 "Couenne: %d cuts (%d row, %d col) for linearization\n",
00638 cs.sizeRowCuts () + cs.sizeColCuts (),
00639 cs.sizeRowCuts (), cs.sizeColCuts ());
00640
00641 fictitiousBound (cs, problem_, true);
00642 firstcall_ = false;
00643 ntotalcuts_ = nrootcuts_ = cs.sizeRowCuts ();
00644
00645 } else {
00646
00647 problem_ -> domain () -> pop ();
00648
00649 ntotalcuts_ += (cs.sizeRowCuts () - nInitCuts);
00650
00651 if (saveOptimum)
00652 realOpt = saveOptimum;
00653 }
00654
00655 septime_ += CoinCpuTime () - now;
00656
00657 if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_CONVEXIFYING)) {
00658
00659 if (cs.sizeColCuts ()) {
00660 jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne col cuts:\n");
00661 for (int i=0; i<cs.sizeColCuts (); i++)
00662 cs.colCutPtr (i) -> print ();
00663 }
00664 }
00665
00666 if (!(info.inTree))
00667 rootTime_ = CoinCpuTime ();
00668 }
00669
00670 }