00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <stdlib.h>
00012
00013 #include "CglCutGenerator.hpp"
00014 #include "CouenneTwoImplied.hpp"
00015 #include "CoinPackedMatrix.hpp"
00016 #include "CouennePrecisions.hpp"
00017 #include "CouenneProblem.hpp"
00018 #include "CouenneExprVar.hpp"
00019
00020 #define MIN_DENOM 1.e-10
00021
00022 namespace Couenne {
00023
00024 enum signum {NEG, POS, DN, UP};
00025
00026
00027
00028
00029
00030 typedef struct {
00031
00032 double alpha;
00033 int indVar;
00034 enum signum sign;
00035
00036
00037 } threshold;
00038
00039
00040
00041
00042 int compthres (register const void *t1,
00043 register const void *t2) {
00044
00045 register double
00046 a1 = (*(threshold **) t1) -> alpha,
00047 a2 = (*(threshold **) t2) -> alpha;
00048
00049 return ((a1 < a2) ? -1 :
00050 (a1 > a2) ? 1 : 0);
00051 }
00052
00053
00054
00055 struct indPosPair {
00056
00057 int index;
00058 int position;
00059 };
00060
00061
00062
00063 int compPair (register const void *p1,
00064 register const void *p2) {
00065
00066 register int
00067 i1 = ((struct indPosPair *) p1) -> index,
00068 i2 = ((struct indPosPair *) p2) -> index;
00069
00070 return ((i1 < i2) ? -1 :
00071 (i1 > i2) ? 1 : 0);
00072 }
00073
00074
00075
00076
00077
00078
00079
00080 int combine (CouenneProblem *p,
00081 int n1,
00082 int n2,
00083 const int *ind1c,
00084 const int *ind2c,
00085 double *sa1,
00086 double *sa2,
00087 const double *a1c,
00088 const double *a2c,
00089 double *clb,
00090 double *cub,
00091 double l1,
00092 double l2,
00093 double u1,
00094 double u2,
00095 bool *isInteger,
00096 int sign) {
00097
00098
00099
00100
00101 int
00102 *ind1 = new int [n1],
00103 *ind2 = new int [n2];
00104
00105 double
00106 *a1 = new double [n1],
00107 *a2 = new double [n2];
00108
00109 CouNumber
00110 *Lb = p -> Lb (),
00111 *Ub = p -> Ub ();
00112
00113 struct indPosPair *pairs = new struct indPosPair [CoinMax (n1,n2)];
00114
00115
00116
00117 for (int i=0; i<n1; i++) {
00118 pairs [i]. index = ind1c [i];
00119 pairs [i]. position = i;
00120 }
00121
00122
00123
00124
00125 qsort (pairs, n1, sizeof (indPosPair), compPair);
00126
00127 for (int i=0; i<n1; i++) {
00128
00129 int rightpos = pairs [i]. position;
00130
00131 ind1 [i] = ind1c [rightpos];
00132 a1 [i] = a1c [rightpos];
00133 }
00134
00136
00137 for (int i=0; i<n2; i++) {
00138 pairs [i]. index = ind2c [i];
00139 pairs [i]. position = i;
00140 }
00141
00142 qsort (pairs, n2, sizeof (indPosPair), compPair);
00143
00144 for (int i=0; i<n2; i++) {
00145
00146 int rightpos = pairs [i]. position;
00147
00148 ind2 [i] = ind2c [rightpos];
00149 a2 [i] = a2c [rightpos];
00150 }
00151
00152 delete [] pairs;
00153
00154
00155
00156
00157 double signA;
00158
00159 if (sign < 0) {
00160
00161 double tmp = u2;
00162 u2 = -l2;
00163 l2 = -tmp;
00164
00165 signA = -1.;
00166
00167 } else signA = 1.;
00168
00169 #ifdef DEBUG
00170 printf ("here are the two ineqs (sign=%d):\n", sign);
00171 printf (" 1: %g <=", l1); for (int i=0; i<n1; i++) printf (" %+g x%d", a1 [i], ind1 [i]);
00172 printf ("<= %g\n 2: %g <=", u1, l2); for (int i=0; i<n2; i++) printf (" %+g x%d", signA * a2 [i], ind2 [i]); printf ("<= %g\n", u2);
00173 #endif
00174
00175 threshold
00176 *alphas = new threshold [n1 + n2],
00177 **inalphas = new threshold* [n1 + n2],
00178 *curalpha = alphas;
00179
00180 int
00181 i1 = 0,
00182 i2 = 0,
00183 cnt = 0,
00184 incnt = 0;
00185
00186 while (i1 < n1 ||
00187 i2 < n2) {
00188
00189 if (i1 < n1 && (i2 == n2 || ind1 [i1] < ind2 [i2])) {
00190
00191 curalpha -> alpha = 1.;
00192 curalpha -> sign = (a1 [i1] < 0.) ? NEG : POS;
00193 curalpha -> indVar = ind1 [i1];
00194
00195 cnt++;
00196 i1++;
00197
00198 } else if (i2 < n2 && (i1 == n1 || ind2 [i2] < ind1 [i1])) {
00199
00200 curalpha -> alpha = 0.;
00201 curalpha -> sign = (signA * a2 [i2] < 0.) ? NEG : POS;
00202 curalpha -> indVar = ind2 [i2];
00203
00204 cnt++;
00205 i2++;
00206
00207 } else {
00208
00209
00210
00211
00212 double a2i2 = signA * a2 [i2];
00213
00214 curalpha -> alpha =
00215 (a1 [i1] == a2i2) ?
00216 -1. :
00217 - a2i2 / (a1 [i1] - a2i2);
00218
00219 if (curalpha -> alpha <= 0. ||
00220 curalpha -> alpha >= 1.)
00221
00222 curalpha -> sign = (a1 [i1] > 0.) ? POS : NEG;
00223
00224 else
00225
00226 curalpha -> sign =
00227 (a1 [i1] == a2i2) ?
00228 ((a1 [i1] < 0.) ? NEG : POS) :
00229 ((a1 [i1] < 0.) ? DN : UP);
00230
00231 curalpha -> indVar = ind1 [i1];
00232
00233 if (curalpha -> alpha > 0. &&
00234 curalpha -> alpha < 1.)
00235 inalphas [incnt++] = curalpha;
00236
00237 cnt++;
00238 i1++;
00239 i2++;
00240 }
00241
00242 #ifdef DEBUG
00243 printf ("(%d,%g,%s) ",
00244 curalpha -> indVar,
00245 curalpha -> alpha,
00246 (curalpha -> sign == NEG) ? "NEG" :
00247 (curalpha -> sign == POS) ? "POS" :
00248 (curalpha -> sign == DN) ? "DN" : "UP");
00249 #endif
00250
00251 curalpha++;
00252 }
00253
00254 #ifdef DEBUG
00255 printf ("\n");
00256 #endif
00257
00258
00259
00260 if (!incnt) {
00261
00262 delete [] inalphas;
00263 delete [] alphas;
00264 delete [] ind1;
00265 delete [] ind2;
00266 delete [] a1;
00267 delete [] a2;
00268
00269 return 0;
00270 }
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281 int ntightened = 0;
00282
00283 if (incnt > 1)
00284 qsort (inalphas, incnt, sizeof (threshold *), compthres);
00285
00286 #ifdef DEBUG
00287 printf ("--sorted: ");
00288 for (int i = 0; i < incnt; i++)
00289 printf ("(%d,%g,%s) ",
00290 inalphas [i] -> indVar,
00291 inalphas [i] -> alpha,
00292 (inalphas [i] -> sign == NEG) ? "NEG" :
00293 (inalphas [i] -> sign == POS) ? "POS" :
00294 (inalphas [i] -> sign == DN) ? "DN" : "UP");
00295 printf ("\n");
00296 #endif
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 double
00321 minSum1 = 0., minSum2 = 0.,
00322 maxSum1 = 0., maxSum2 = 0.;
00323
00324 int
00325 mInfs = 0,
00326 pInfs = 0;
00327
00328
00329
00330
00331
00332
00333
00334 for (int i=0; i<n1; i++) {
00335
00336 int indVar1 = ind1 [i];
00337
00338 double
00339 a1i = a1 [i],
00340 a2i = sa2 [indVar1],
00341 clb1 = clb [indVar1],
00342 cub1 = cub [indVar1];
00343
00344
00345
00346
00347 if (a2i == 0.)
00348 a2i = a1i;
00349
00350 #ifdef DEBUG
00351 else if (fabs (a2i) < 1e-10)
00352 printf ("\nNumerics in a2i\n\n");
00353 #endif
00354
00355 if (a2i < 0.) {
00356
00357 bool zeroA = (fabs (sa2 [indVar1]) == 0.);
00358
00359 if (cub1 > COUENNE_INFINITY/10) {if (zeroA) mInfs ++;} else minSum1 += a1i * cub1;
00360 if (clb1 < - COUENNE_INFINITY/10) {if (zeroA) pInfs ++;} else maxSum1 += a1i * clb1;
00361
00362 } else if (a2i > 0.) {
00363
00364 bool zeroA = (fabs (sa2 [indVar1]) == 0.);
00365
00366 if (clb1 < - COUENNE_INFINITY/10) {if (zeroA) mInfs ++;} else minSum1 += a1i * clb1;
00367 if (cub1 > COUENNE_INFINITY/10) {if (zeroA) pInfs ++;} else maxSum1 += a1i * cub1;
00368 }
00369 }
00370
00371
00372
00373 for (int i=0; i<n2; i++) {
00374
00375 int indVar2 = ind2 [i];
00376
00377 double
00378 a2i = a2 [i] * signA,
00379 clb2 = clb [indVar2],
00380 cub2 = cub [indVar2];
00381
00382 if (a2i < 0.) {
00383
00384 if (cub2 > COUENNE_INFINITY/10) mInfs ++; else minSum2 += a2i * cub2;
00385 if (clb2 < - COUENNE_INFINITY/10) pInfs ++; else maxSum2 += a2i * clb2;
00386
00387 } else {
00388
00389 if (clb2 < - COUENNE_INFINITY/10) mInfs ++; else minSum2 += a2i * clb2;
00390 if (cub2 > COUENNE_INFINITY/10) pInfs ++; else maxSum2 += a2i * cub2;
00391 }
00392 }
00393
00394
00395
00396
00397
00398
00399
00400 #ifdef DEBUG
00401 printf (" at alpha=zero, m[in,ax]sum[12] = (1:(%g,%g), 2:(%g,%g)), %d mInf, %d pInf\n",
00402 minSum1, maxSum1,
00403 minSum2, maxSum2,
00404 mInfs, pInfs);
00405 #endif
00406
00407
00408
00409
00410
00411 std::vector <std::pair <int, double> >
00412 newLB,
00413 newUB;
00414
00415 delete [] a1;
00416 delete [] a2;
00417
00418 #ifdef DEBUG
00419 printf (" ");
00420 for (int j=0; j < n1; j++) printf ("x%d [%g,%g] ", ind1 [j], clb [ind1 [j]], cub [ind1 [j]]); printf ("--- ");
00421 for (int j=0; j < n2; j++) printf ("x%d [%g,%g] ", ind2 [j], clb [ind2 [j]], cub [ind2 [j]]); printf ("\n");
00422 #endif
00423
00424 for (int i = 0; i < incnt; i++) {
00425
00426 threshold *inalpha = inalphas [i];
00427
00428 #ifdef DEBUG
00429 printf (" looking at %d: (%d,%g,%s)\n", i,
00430 inalphas [i] -> indVar,
00431 inalphas [i] -> alpha,
00432 (inalphas [i] -> sign == NEG) ? "NEG" :
00433 (inalphas [i] -> sign == POS) ? "POS" :
00434 (inalphas [i] -> sign == DN) ? "DN" : "UP");
00435
00436 fflush (stdout);
00437 #endif
00438
00439
00440
00441
00442
00443
00444
00445 double alpha = inalpha -> alpha;
00446
00447
00448
00449
00450 while (true) {
00451
00452 inalpha = inalphas [i];
00453
00454 int
00455 signalpha = inalpha -> sign,
00456 indVar = inalpha -> indVar;
00457
00458 double
00459 clbi = clb [indVar],
00460 cubi = cub [indVar];
00461
00462
00463
00464
00465 if (fabs (clbi) > 0.) {
00466
00467 if (clbi < - COUENNE_INFINITY/10) {
00468
00469 if (signalpha == DN) {pInfs++; mInfs--;}
00470 else {pInfs--; mInfs++;}
00471
00472 } else {
00473
00474 double
00475 sa1ic = sa1 [indVar] * clbi,
00476 sa2ic = sa2 [indVar] * clbi;
00477
00478 if (signalpha == DN) {
00479
00480 minSum1 -= sa1ic; minSum2 -= sa2ic;
00481 maxSum1 += sa1ic; maxSum2 += sa2ic;
00482
00483 } else {
00484
00485 minSum1 += sa1ic; minSum2 += sa2ic;
00486 maxSum1 -= sa1ic; maxSum2 -= sa2ic;
00487 }
00488 }
00489 }
00490
00491
00492
00493 if (fabs (cubi) > 0.) {
00494
00495 if (cubi > COUENNE_INFINITY/10) {
00496
00497 if (signalpha == DN) {pInfs--; mInfs++;}
00498 else {pInfs++; mInfs--;}
00499
00500 } else {
00501
00502 double
00503 sa1ic = sa1 [indVar] * cubi,
00504 sa2ic = sa2 [indVar] * cubi;
00505
00506 if (signalpha == DN) {
00507
00508 minSum1 += sa1ic; minSum2 += sa2ic;
00509 maxSum1 -= sa1ic; maxSum2 -= sa2ic;
00510
00511 } else {
00512
00513 minSum1 -= sa1ic; minSum2 -= sa2ic;
00514 maxSum1 += sa1ic; maxSum2 += sa2ic;
00515 }
00516 }
00517 }
00518
00519 if ((i < incnt-1) &&
00520 (inalphas [i+1] -> alpha == alpha)) ++i;
00521 else break;
00522 }
00523
00524 #ifdef DEBUG
00525 printf (" at alpha=%g, m[in,ax]sum[12] = (1:(%g,%g), 2:(%g,%g)), %d mInf, %d pInf\n",
00526 inalphas [i] -> alpha,
00527 minSum1, maxSum1,
00528 minSum2, maxSum2,
00529 mInfs, pInfs);
00530 #endif
00531
00532 for (int j=0; j<n1+n2; j++) {
00533
00534
00535
00536
00537
00538
00539
00540 if (j >= n1 && sa1 [ind2 [j - n1]] != 0.)
00541 continue;
00542
00543 int
00544 indVar = ((j < n1) ? ind1 [j] : ind2 [j - n1]),
00545 tickMin = 0,
00546 tickMax = 0;
00547
00548 double
00549 sa1i = sa1 [indVar],
00550 sa2i = sa2 [indVar],
00551 ci = alpha * sa1i + (1. - alpha) * sa2i;
00552
00553
00554
00555
00556 if (fabs (ci) < MIN_DENOM)
00557 continue;
00558
00559 double
00560 clbi = clb [indVar],
00561 cubi = cub [indVar],
00562
00563 newL = - COUENNE_INFINITY,
00564 newU = COUENNE_INFINITY,
00565
00566 subMin1 = 0., subMin2 = 0.,
00567 subMax1 = 0., subMax2 = 0.;
00568
00569 #ifdef DEBUG
00570 printf (" now tightening x%d\n denominator: %g * %g + %g * %g = %g\n",
00571 (j < n1) ? ind1 [j] : ind2 [j-n1],
00572 alpha, sa1i, 1-alpha, sa2i, ci);
00573 #endif
00574
00575
00576
00577
00578
00579 if (clbi <= - COUENNE_INFINITY/10) {
00580
00581
00582
00583
00584
00585 if (ci > 0.) tickMin++;
00586 else tickMax++;
00587
00588 } else {
00589
00590 if (ci > 0.) {
00591
00592 subMin1 = sa1i * clbi;
00593 subMin2 = sa2i * clbi;
00594
00595 } else {
00596
00597 subMax1 = sa1i * clbi;
00598 subMax2 = sa2i * clbi;
00599 }
00600 }
00601
00602 if (cubi >= COUENNE_INFINITY/10) {
00603
00604
00605
00606
00607
00608 if (ci > 0.) tickMax++;
00609 else tickMin++;
00610
00611 } else {
00612
00613 if (ci > 0.) {
00614
00615 subMax1 = sa1i * cubi;
00616 subMax2 = sa2i * cubi;
00617
00618 } else {
00619
00620 subMin1 = sa1i * cubi;
00621 subMin2 = sa2i * cubi;
00622 }
00623 }
00624
00625 if (fabs (ci) < MIN_DENOM) {
00626
00627 } else {
00628
00629 if ((l1 > -COUENNE_INFINITY/10) &&
00630 (l2 > -COUENNE_INFINITY/10) &&
00631 (pInfs == tickMax)) {
00632
00633 newL = ((l1 - l2 - (maxSum1 - maxSum2) + (subMax1 - subMax2)) * alpha + l2 - maxSum2 + subMax2) / ci;
00634
00635 #ifdef DEBUG
00636 printf (" \
00637 attempting newL = ((l1 - l2 - (maxSum1 - maxSum2) + (subMax1 - subMax2)) * alpha + l2 - maxSum2 + subMax2) / ci\n \
00638 ((%g - %g - (%g - %g) + (%g - %g)) * %g + %g - %g + %g) / %g = %g\n",
00639 l1, l2, maxSum1, maxSum2, subMax1, subMax2, alpha, l2, maxSum2, subMax2, ci, newL);
00640 #endif
00641 }
00642
00643 if ((u1 < COUENNE_INFINITY/10) &&
00644 (u2 < COUENNE_INFINITY/10) &&
00645 (mInfs == tickMin)) {
00646
00647 newU = ((u1 - u2 - (minSum1 - minSum2) + (subMin1 - subMin2)) * alpha + u2 - minSum2 + subMin2) / ci;
00648
00649 #ifdef DEBUG
00650 printf (" \
00651 attempting newU = ((u1 - u2 - (minSum1 - minSum2) + (subMin1 - subMin2)) * alpha + u2 - minSum2 + subMin2) / ci\n \
00652 ((%g - %g - (%g - %g) + (%g - %g)) * %g + %g - %g + %g) / %g = %g\n",
00653 u1, u2, minSum1, minSum2, subMin1, subMin2, alpha, u2, minSum2, subMin2, ci, newU);
00654 #endif
00655 }
00656
00657 if (ci < 0.) {
00658
00659 #ifdef DEBUG
00660 printf (" swap'em: %g, %g\n", newL, newU);
00661 #endif
00662
00663 register double tmp = newL <= - COUENNE_INFINITY / 10 ? COUENNE_INFINITY : newL;
00664 newL = newU >= COUENNE_INFINITY / 10 ? - COUENNE_INFINITY : newU;
00665 newU = tmp;
00666 }
00667
00668 #ifdef DEBUG
00669 printf (" final: %g, %g\n", newL, newU);
00670 #endif
00671 }
00672
00673 #ifdef DEBUG
00674 printf (" bounds for x_%d: [%g,%g] vs [%g,%g]\n", indVar, newL, newU, clb [indVar], cub [indVar]);
00675 #endif
00676
00677 if ((newL > cubi + COUENNE_EPS) ||
00678 (newU < clbi - COUENNE_EPS)) {
00679
00680 delete [] inalphas;
00681 delete [] alphas;
00682 delete [] ind1;
00683 delete [] ind2;
00684
00685 #ifdef DEBUG
00686 printf ("infeasible\n");
00687 #endif
00688 return -1;
00689 }
00690
00691 #ifdef DEBUG
00692 if (p -> bestSol () &&
00693 ((p -> bestSol () [indVar] > newU + COUENNE_EPS) ||
00694 (p -> bestSol () [indVar] < newL - COUENNE_EPS)) &&
00695 (p -> bestSol () [indVar] >= clbi) &&
00696 (p -> bestSol () [indVar] <= cubi))
00697
00698 printf ("optimum violated: %g not in [%g,%g]\n", p -> bestSol () [indVar], newL, newU);
00699 #endif
00700
00701 double
00702 &lbi = Lb [indVar],
00703 &ubi = Ub [indVar];
00704
00705 if ((newL > lbi + COUENNE_EPS) && (newL > -COUENNE_INFINITY / 10)) {
00706
00707 ntightened++;
00708 lbi = newL;
00709 newLB. push_back (std::pair <int, double> (indVar, newL));
00710
00711 #ifdef DEBUG
00712 printf (" new lower bound: x%d >= %g [%g]\n", indVar, newL, clb [indVar]);
00713 #endif
00714 }
00715
00716 if ((newU < ubi - COUENNE_EPS) && (newU < COUENNE_INFINITY / 10)) {
00717
00718 ntightened++;
00719 ubi = newU;
00720 newUB. push_back (std::pair <int, double> (indVar, newU));
00721
00722 #ifdef DEBUG
00723 printf (" new upper bound: x%d <= %g [%g]\n", indVar, newU, cub [indVar]);
00724 #endif
00725 }
00726 }
00727 }
00728
00729 #ifdef DEBUG
00730 printf ("===================================\n");
00731 #endif
00732
00733 for (std::vector <std::pair <int, double> >::iterator i = newLB. begin (); i != newLB. end (); ++i) {
00734
00735 if (isInteger [i -> first]) i -> second = ceil (i -> second - COUENNE_EPS);
00736
00737 if (i -> second > clb [i -> first])
00738 Lb [i -> first] = clb [i -> first] = i -> second;
00739 }
00740
00741 for (std::vector <std::pair <int, double> >::iterator i = newUB. begin (); i != newUB. end (); ++i) {
00742
00743 if (isInteger [i -> first]) i -> second = floor (i -> second + COUENNE_EPS);
00744
00745 if (i -> second < cub [i -> first])
00746 Ub [i -> first] = cub [i -> first] = i -> second;
00747 }
00748
00749 delete [] inalphas;
00750 delete [] alphas;
00751
00752 delete [] ind1;
00753 delete [] ind2;
00754
00755 return ntightened;
00756 }
00757
00758 }