00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <set>
00012 #include <stdlib.h>
00013
00014 #include "BonCbc.hpp"
00015 #include "BonBabInfos.hpp"
00016 #include "CoinHelperFunctions.hpp"
00017 #include "CoinPackedMatrix.hpp"
00018 #include "CglCutGenerator.hpp"
00019 #include "CoinTime.hpp"
00020
00021 #include "CouenneTypes.hpp"
00022 #include "CouenneProblemElem.hpp"
00023 #include "CouenneTwoImplied.hpp"
00024 #include "CouenneExprVar.hpp"
00025 #include "CouennePrecisions.hpp"
00026 #include "CouenneProblem.hpp"
00027 #include "CouenneInfeasCut.hpp"
00028 #include "CouenneJournalist.hpp"
00029
00030 using namespace Ipopt;
00031
00032
00033 namespace Couenne {
00034
00036 void updateBranchInfo (const OsiSolverInterface &si, CouenneProblem *p,
00037 t_chg_bounds *chg, const CglTreeInfo &info);
00038
00039
00040
00041 int combine (CouenneProblem *p,
00042 int n1, int n2,
00043 const int *ind1,
00044 const int *ind2,
00045 double *sa1,
00046 double *sa2,
00047 const double *a1,
00048 const double *a2,
00049 double *clb,
00050 double *cub,
00051 double l1,
00052 double l2,
00053 double u1,
00054 double u2,
00055 bool *isInteger,
00056 int sign);
00057
00059 void CouenneTwoImplied::generateCuts (const OsiSolverInterface &si,
00060 OsiCuts &cs,
00061 const CglTreeInfo info)
00062 #if CGL_VERSION_MAJOR == 0 && CGL_VERSION_MINOR <= 57
00063 const
00064 #endif
00065 {
00066
00067
00068
00069
00070
00071 if (isWiped (cs))
00072 return;
00073
00074 double now = CoinCpuTime ();
00075
00076
00077
00078 if ((depthStopSeparate_ >= 0 &&
00079 info.level > depthStopSeparate_)
00080 ||
00081 (depthLevelling_ >= 0 &&
00082 info.level >= depthLevelling_ &&
00083 CoinDrand48 () > 1. / (2. + info.level - depthLevelling_)))
00084 return;
00085
00086 if (info.level <= 0)
00087 jnlst_ -> Printf (J_ERROR, J_COUENNE, "TwoImpl-BT: "); fflush (stdout);
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 problem_ -> domain () -> push (&si, &cs);
00099
00100 static int nBadColMatWarnings = 0;
00101
00102 std::set <std::pair <int, int> > pairs;
00103
00110
00113
00114 #define USE_ROW
00115
00116 #ifdef USE_ROW
00117
00118 const CoinPackedMatrix *mat = si. getMatrixByRow ();
00119
00120 int
00121 m = mat -> getMajorDim (),
00122 n = mat -> getMinorDim ();
00123
00124 #else
00125
00126 const CoinPackedMatrix *mat = si. getMatrixByCol ();
00127
00128 int
00129 n = mat -> getMajorDim (),
00130 m = mat -> getMinorDim ();
00131
00132 #endif
00133
00134 const double
00135 *rlb = si.getRowLower (),
00136 *rub = si.getRowUpper ();
00137
00138 #ifdef USE_ROW
00139
00141
00142 int
00143 nnz = mat -> getNumElements (),
00144 nnzC = 0,
00145 *sta = new int [n+1],
00146 nCuts = cs.sizeRowCuts ();
00147
00148
00149
00150 for (int i=0, ii = cs. sizeRowCuts (); ii--; i++) {
00151
00152 const OsiRowCut *cut = cs. rowCutPtr (i);
00153 const CoinPackedVector &rowCoe = cut -> row ();
00154
00155 nnzC += rowCoe.getNumElements ();
00156 }
00157
00158 int *ind = new int [nnz + nnzC];
00159 double *A = new double [nnz + nnzC];
00160
00162 {
00163 const double
00164 *rA = mat -> getElements ();
00165
00166 const int
00167 *rInd = mat -> getIndices (),
00168 *rSta = mat -> getVectorStarts ();
00169
00170
00171
00172 CoinZeroN (sta, n+1);
00173
00175
00176
00177
00178 for (int i=nnz; i--;)
00179 ++ (sta [1 + *rInd++]);
00180
00181
00182
00183 for (int i=0, ii = cs. sizeRowCuts (); ii--; i++) {
00184
00185 const OsiRowCut *cut = cs. rowCutPtr (i);
00186 const CoinPackedVector &rowCoe = cut -> row ();
00187 const int *indices = rowCoe.getIndices ();
00188 int nnz = rowCoe.getNumElements ();
00189
00190
00191 for (int i=nnz; i--;)
00192 ++ (sta [1 + *indices++]);
00193 }
00194
00195 rInd -= nnz;
00196
00198
00199
00200
00201 for (int i=1; i<=n; i++)
00202 sta [i] += sta [i-1];
00203
00204
00205
00206
00207 for (int i=0; i<m; i++) {
00208
00209
00210
00211 int rowStart = rSta [i];
00212
00213 for (int j = rowStart, jj = rSta [i+1] - rowStart; jj--; j++) {
00214
00215 int &curSta = sta [rInd [j]];
00216
00217 ind [curSta] = i;
00218 A [curSta++] = rA [j];
00219 }
00220
00221
00222 }
00223
00224
00225
00226 for (int i=0, ii = cs. sizeRowCuts (); ii--; i++) {
00227
00228 const OsiRowCut *cut = cs. rowCutPtr (i);
00229
00230
00231 const CoinPackedVector &rowCoe = cut -> row ();
00232 const int *indices = rowCoe.getIndices ();
00233 const double *elements = rowCoe.getElements ();
00234 int nnz = rowCoe.getNumElements ();
00235
00236
00237 for (int j=nnz; j--;) {
00238
00239
00240
00241 int &curSta = sta [*indices++];
00242
00243 ind [curSta] = m+i;
00244 A [curSta++] = *elements++;
00245 }
00246
00247 }
00248
00249 for (int i=n; --i;)
00250 sta [i] = sta [i-1];
00251
00252 sta [0] = 0;
00253
00254
00255 }
00256
00257 #else
00258
00259 const double
00260 *A = mat -> getElements ();
00261
00262 const int
00263 *ind = mat -> getIndices (),
00264 *sta = mat -> getVectorStarts ();
00265
00266 #endif
00267
00270
00271 bool *isInteger = new bool [n];
00272 for (int i=0, ii=n; ii--; i++)
00273 *isInteger++ = problem_ -> Var (i) -> isInteger ();
00274 isInteger -= n;
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298 double
00299 *sa1 = new double [n],
00300 *sa2 = new double [n];
00301
00302 CoinFillN (sa1, n, 0.);
00303 CoinFillN (sa2, n, 0.);
00304
00305 for (int i=0; i<n; i++, sta++) {
00306
00307
00308
00309 for (int jj = *(sta+1) - *sta, j = *sta; jj--; j++)
00310 for (int kk = jj, k = j+1; kk--; k++) {
00311
00312 if (CoinCpuTime () > problem_ -> getMaxCpuTime ())
00313 break;
00314
00315 register int
00316 indj = ind [j],
00317 indk = ind [k];
00318
00319
00320
00321
00322 if ((indj >= m + nCuts) || (indj < 0) ||
00323 (indk >= m + nCuts) || (indk < 0)) {
00324
00325 if (nBadColMatWarnings++ < 1)
00326
00327 printf ("\
00328 Couenne: warning, matrix by row has nonsense indices.\n\
00329 This separator will now return without (column) cuts.\n\
00330 NOTE: further such inconsistencies won't be reported.\n");
00331
00332 delete [] sa1;
00333 delete [] sa2;
00334
00335 delete [] sta;
00336 delete [] ind;
00337 delete [] A;
00338
00339 delete [] isInteger;
00340
00341 problem_ -> domain () -> pop ();
00342
00343 totalTime_ += CoinCpuTime () - now;
00344 totalInitTime_ += CoinCpuTime () - now;
00345
00346 return;
00347 }
00348
00349 double
00350 prod = A [j] * A [k],
00351 rlbj, rubj, rlbk, rubk;
00352
00353 if (indj>=m) {
00354 OsiRowCut *cut = cs.rowCutPtr (indj-m);
00355 rlbj = cut -> lb ();
00356 rubj = cut -> ub ();
00357 } else {
00358 rlbj = rlb [indj];
00359 rubj = rub [indj];
00360 }
00361
00362 if (indk>=m) {
00363 OsiRowCut *cut = cs.rowCutPtr (indk-m);
00364 rlbk = cut -> lb ();
00365 rubk = cut -> ub ();
00366 } else {
00367 rlbk = rlb [indk];
00368 rubk = rub [indk];
00369 }
00370
00371 if (prod > 0.) {
00372
00373
00374
00375
00376 if (!(((rlbj > -COUENNE_INFINITY/10) && (rubk < COUENNE_INFINITY/10)) ||
00377 ((rubj < COUENNE_INFINITY/10) && (rlbk > -COUENNE_INFINITY/10))))
00378
00379 continue;
00380
00381 } else
00382
00383 if ((prod < 0.) &&
00384
00385 !(((rlbj > -COUENNE_INFINITY/10) && (rlbk > -COUENNE_INFINITY/10)) ||
00386 ((rubj < COUENNE_INFINITY/10) && (rubk < COUENNE_INFINITY/10))))
00387
00388 continue;
00389
00390 pairs.insert (std::pair <int, int> (indj, indk));
00391 }
00392 }
00393
00394
00395
00396 #ifdef USE_ROW
00397 const CoinPackedMatrix *rowA = mat;
00398 #else
00399 const CoinPackedMatrix *rowA = si. getMatrixByRow ();
00400 #endif
00401
00402 const double
00403 *rA = rowA -> getElements ();
00404
00405
00406
00407 double
00408 *clb = CoinCopyOfArray (problem_ -> Lb (), n),
00409 *cub = CoinCopyOfArray (problem_ -> Ub (), n),
00410 *oldLB = CoinCopyOfArray (problem_ -> Lb (), n),
00411 *oldUB = CoinCopyOfArray (problem_ -> Ub (), n);
00412
00413 const int
00414 *rInd = rowA -> getIndices (),
00415 *rSta = rowA -> getVectorStarts ();
00416
00417 int
00418 ntightened = 0,
00419 ntrials = 0,
00420 nCurTightened;
00421
00422
00423
00424 Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (si.getAuxiliaryInfo ());
00425
00426 int result = 0;
00427
00428
00429
00430 t_chg_bounds *chg_bds = new t_chg_bounds [n];
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440 updateBranchInfo (si, problem_, chg_bds, info);
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 do {
00454
00455 nCurTightened = 0;
00456
00457
00458
00459
00460 for (std::set <std::pair <int, int> >:: iterator p = pairs.begin (); p != pairs.end (); ++p) {
00461
00462 if (CoinCpuTime () > problem_ -> getMaxCpuTime ())
00463 break;
00464
00465
00466
00467 int
00468 h = p -> first,
00469 k = p -> second,
00470 n1, n2;
00471
00472 double l1, u1, l2, u2;
00473
00474 const double *a1, *a2;
00475
00476 const int *ind1, *ind2;
00477
00478
00479
00480
00481 if (h>=m) {
00482
00483 const OsiRowCut *cut = cs.rowCutPtr (h-m);
00484 const CoinPackedVector &rowCoe = cut -> row ();
00485
00486 l1 = cut -> lb ();
00487 u1 = cut -> ub ();
00488 n1 = rowCoe. getNumElements ();
00489 ind1 = rowCoe. getIndices ();
00490 a1 = rowCoe. getElements ();
00491
00492 } else {
00493
00494 l1 = rlb [h];
00495 u1 = rub [h];
00496 n1 = rSta [h+1] - rSta [h];
00497 ind1 = rInd + rSta [h];
00498 a1 = rA + rSta [h];
00499 }
00500
00502
00503 if (k>=m) {
00504
00505 OsiRowCut *cut = cs.rowCutPtr (k-m);
00506 const CoinPackedVector &rowCoe = cut -> row ();
00507
00508 l2 = cut -> lb ();
00509 u2 = cut -> ub ();
00510 n2 = rowCoe. getNumElements ();
00511 ind2 = rowCoe. getIndices ();
00512 a2 = rowCoe. getElements ();
00513
00514 } else {
00515
00516 l2 = rlb [k];
00517 u2 = rub [k];
00518 n2 = rSta [k+1] - rSta [k];
00519 ind2 = rInd + rSta [k];
00520 a2 = rA + rSta [k];
00521 }
00522
00523
00524
00525
00526
00527 if ((h < m) && (k < m) && !firstCall_) {
00528
00529 bool new_bound = false;
00530
00531 for (int i=n1; i--;) {
00532
00533 t_chg_bounds &chg = chg_bds [ind1 [i]];
00534
00535 if ((chg. lower () != t_chg_bounds::UNCHANGED) ||
00536 (chg. upper () != t_chg_bounds::UNCHANGED))
00537
00538 new_bound = true;
00539 }
00540
00541 if (!new_bound)
00542 for (int i=n2; i--;) {
00543
00544 t_chg_bounds &chg = chg_bds [ind2 [i]];
00545
00546 if ((chg. lower () != t_chg_bounds::UNCHANGED) ||
00547 (chg. upper () != t_chg_bounds::UNCHANGED))
00548
00549 new_bound = true;
00550 }
00551
00552 if (!new_bound) continue;
00553 }
00554
00555
00556
00557 for (int i=n1; i--;) sa1 [ind1 [i]] = a1 [i];
00558 for (int i=n2; i--;) sa2 [ind2 [i]] = a2 [i];
00559
00560
00561
00562
00563
00564
00565
00566 if ((u1 < COUENNE_INFINITY && u2 < COUENNE_INFINITY) ||
00567 (l1 > - COUENNE_INFINITY && l2 > - COUENNE_INFINITY))
00568
00569 result = combine (problem_, n1, n2, ind1, ind2, sa1, sa2, a1, a2, clb, cub, l1, l2, u1, u2, isInteger, 1);
00570
00571 if (result < 0)
00572 break;
00573
00574 nCurTightened += result;
00575 result = 0;
00576
00577
00578 for (int i=n2; i--;)
00579 sa2 [ind2 [i]] = - a2 [i];
00580
00581 if ((u1 < COUENNE_INFINITY && l2 > - COUENNE_INFINITY) ||
00582 (l1 > - COUENNE_INFINITY && u2 < COUENNE_INFINITY))
00583
00584
00585 result = combine (problem_, n1, n2, ind1, ind2, sa1, sa2, a1, a2, clb, cub, l1, l2, u1, u2, isInteger, -1);
00586
00587 if (result < 0)
00588 break;
00589
00590 nCurTightened += result;
00591
00592
00593
00594 for (int i=n1; i--;) sa1 [ind1 [i]] = 0.;
00595 for (int i=n2; i--;) sa2 [ind2 [i]] = 0.;
00596 }
00597
00598 if (result < 0)
00599 break;
00600
00601 int objInd = problem_ -> Obj (0) -> Body () -> Index ();
00602
00603 if (nCurTightened &&
00604 (objInd >= 0) &&
00605 babInfo &&
00606 babInfo -> babPtr ()) {
00607
00608 #ifdef DEBUG
00609 printf ("FBBT\n");
00610 #endif
00611
00612 CouNumber
00613 UB = babInfo -> babPtr () -> model (). getObjValue(),
00614 LB = babInfo -> babPtr () -> model (). getBestPossibleObjValue (),
00615 primal0 = problem_ -> Ub (objInd),
00616 dual0 = problem_ -> Lb (objInd);
00617
00618
00619
00620 if ((UB < COUENNE_INFINITY) &&
00621 (UB < primal0 - COUENNE_EPS)) {
00622
00623 problem_ -> Ub (objInd) = UB;
00624 chg_bds [objInd].setUpper (t_chg_bounds::CHANGED);
00625 }
00626
00627 if ((LB > - COUENNE_INFINITY) &&
00628 (LB > dual0 + COUENNE_EPS)) {
00629 problem_ -> Lb (objInd) = LB;
00630 chg_bds [objInd].setLower (t_chg_bounds::CHANGED);
00631 }
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657 if (!(problem_ -> btCore (chg_bds))) {
00658
00659
00660
00661 result = -1;
00662 break;
00663
00664 } else {
00665
00666
00667
00668 for (int i=0; i<n; i++) {
00669
00670 if (problem_ -> Lb (i) > clb [i] + COUENNE_EPS) clb [i] = problem_ -> Lb (i);
00671 if (problem_ -> Ub (i) < cub [i] - COUENNE_EPS) cub [i] = problem_ -> Ub (i);
00672 }
00673 }
00674 }
00675
00676 ntightened += nCurTightened;
00677
00678 } while (++ntrials < nMaxTrials_ && nCurTightened);
00679
00680 totalInitTime_ += CoinCpuTime () - now;
00681
00682
00683
00684 if (result >= 0 && ntightened) {
00685
00686
00687
00688 int
00689 *indLB = new int [n],
00690 *indUB = new int [n],
00691 ntightenedL = 0,
00692 ntightenedU = 0;
00693
00694 double
00695 *valLB = new double [n],
00696 *valUB = new double [n];
00697
00698 for (int i=0; i<n; i++) {
00699
00700 if (clb [i] > oldLB [i]) {
00701
00702 if (problem_ -> bestSol () &&
00703 problem_ -> bestSol () [i] > oldLB [i] &&
00704 problem_ -> bestSol () [i] < clb [i] - 1e-12) {
00705
00706 jnlst_ -> Printf (J_ERROR, J_COUENNE,
00707 "Warning, twoImplBounds new LB cuts optimal solution: LB x_%d = %g --> %g, opt %g (diff: %g)\n",
00708 i, oldLB [i], clb [i], problem_ -> bestSol () [i], clb [i] - problem_ -> bestSol () [i]);
00709
00710 }
00711
00712
00713 indLB [ntightenedL] = i;
00714 valLB [ntightenedL++] = clb [i];
00715 }
00716
00717 if (cub [i] < oldUB [i]) {
00718
00719 if (problem_ -> bestSol () &&
00720 problem_ -> bestSol () [i] < oldUB [i] &&
00721 problem_ -> bestSol () [i] > cub [i] + 1e-12) {
00722
00723 jnlst_ -> Printf (J_ERROR, J_COUENNE,
00724 "Warning, twoImplBounds new UB cuts optimal solution: UB x_%d = %g --> %g, opt %g (diff: %g)\n",
00725 i, oldUB [i], cub [i], problem_ -> bestSol () [i], problem_ -> bestSol () [i] - cub [i]);
00726
00727 }
00728
00729
00730 indUB [ntightenedU] = i;
00731 valUB [ntightenedU++] = cub [i];
00732 }
00733 }
00734
00735 if (ntightenedL || ntightenedU) {
00736
00737 OsiColCut newBound;
00738
00739 newBound.setLbs (ntightenedL, indLB, valLB);
00740 newBound.setUbs (ntightenedU, indUB, valUB);
00741
00742
00743
00744 cs.insert (newBound);
00745 }
00746
00747 delete [] indLB;
00748 delete [] indUB;
00749 delete [] valLB;
00750 delete [] valUB;
00751
00752 ntightened = ntightenedL + ntightenedU;
00753 }
00754
00755 else
00756
00757 if (result < 0)
00758 WipeMakeInfeas (cs);
00759
00760 delete [] clb;
00761 delete [] cub;
00762 delete [] oldLB;
00763 delete [] oldUB;
00764 delete [] sa1;
00765 delete [] sa2;
00766 delete [] chg_bds;
00767
00768 delete [] isInteger;
00769
00770 problem_ -> domain () -> pop ();
00771
00772 #ifdef USE_ROW
00773 delete [] A;
00774 delete [] ind;
00775 delete [] (sta-n);
00776 #endif
00777
00778 if (firstCall_)
00779 firstCall_ = false;
00780
00781 totalTime_ += CoinCpuTime () - now;
00782
00783 if (info.level <= 0)
00784 jnlst_ -> Printf (J_ERROR, J_COUENNE, "%d improved bounds\n", ntightened);
00785 }
00786 }