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