00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CoinHelperFunctions.hpp"
00012
00013 #include "CouenneProblem.hpp"
00014 #include "CouenneProblemElem.hpp"
00015 #include "CouenneExprVar.hpp"
00016
00017 #include "BonBabInfos.hpp"
00018
00019 using namespace Couenne;
00020
00021 #define MAX_ITER 3 // max # fake tightening (inner) iterations
00022 #define AGGR_MUL 2 // the larger, the more conservative. Must be > 0
00023 #define AGGR_DIV 2 // the smaller, the more conservative. Must be > 1
00024
00025
00026 const CouNumber phi = 0.5 * (1. + sqrt (5.));
00027
00028
00029 CouNumber fictBounds (char direction,
00030 CouNumber x,
00031 CouNumber lb,
00032 CouNumber ub) {
00033
00034 #define LARGE_BOUND 1e10
00035
00036 if (lb < -LARGE_BOUND) {
00037 if (ub > LARGE_BOUND) {
00038
00039 return (!direction ? -sqrt (-lb) : sqrt (ub));
00040
00041
00042
00043
00044 } else {
00045
00046 if (!direction)
00047 return -sqrt (-lb);
00048
00049 if (x < -COUENNE_EPS) return (CoinMin (0., (x+ub)/2));
00050 else if (x > COUENNE_EPS) return ((x + (ub-x)/AGGR_DIV));
00051 else return ((ub/AGGR_DIV));
00052
00053
00054
00055
00056 }
00057 }
00058 else {
00059 if (ub > LARGE_BOUND) {
00060
00061 if (direction)
00062 return sqrt (ub);
00063
00064 if (x < -COUENNE_EPS) return ((x - (x-lb) / AGGR_DIV));
00065 else if (x > COUENNE_EPS) return (CoinMax (0.,(x+lb)/2));
00066 else return (lb/AGGR_DIV);
00067
00068
00069
00070
00071
00072 } else
00073 return (direction ?
00074 (x + (ub-x) / AGGR_DIV) :
00075 (x - (x-lb) / AGGR_DIV));
00076 }
00077 }
00078
00079
00080
00081
00082
00083
00084
00085 int CouenneProblem::
00086 fake_tighten (char direction,
00087 int index,
00088 const double *X,
00089 CouNumber *olb,
00090 CouNumber *oub,
00091 t_chg_bounds *chg_bds,
00092 t_chg_bounds *f_chg) const {
00093 int
00094 ncols = nVars (),
00095 objind = Obj (0) -> Body () -> Index ();
00096
00097 assert (objind >= 0);
00098
00099 bool
00100 tightened = false,
00101 intvar = variables_ [index] -> isInteger ();
00102
00103 CouNumber
00104 xcur = X [index],
00105 inner = xcur,
00106 outer = (direction ? oub : olb) [index],
00107 fb = fictBounds (direction, xcur, Lb (index), Ub (index));
00108
00109
00110
00111
00112
00113 jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING,
00114 " x_%d. x = %10g, lb = %g, cutoff = %g-----------------\n",
00115 index,xcur,Lb (objind),getCutOff());
00116
00117
00118
00119
00120
00121 for (int iter = 0; iter < MAX_ITER; iter++) {
00122
00123 if (intvar) {
00124
00125 if (!direction) {inner = floor (inner + COUENNE_EPS); outer = ceil (outer - COUENNE_EPS);}
00126 else {inner = ceil (inner - COUENNE_EPS); outer = floor (outer + COUENNE_EPS);}
00127
00128 if ( (direction && (inner > outer + .5)) ||
00129 (!direction && (inner < outer - .5))) {
00130
00131
00132
00133
00134 if (direction) {oub[index] = Ub (index) = fb; chg_bds[index].setUpper(t_chg_bounds::CHANGED);}
00135 else {olb[index] = Lb (index) = fb; chg_bds[index].setLower(t_chg_bounds::CHANGED);}
00136
00137 tightened = true;
00138
00139 if (!(btCore (f_chg)))
00140 return -1;
00141
00142 CoinCopyN (Lb (), ncols, olb);
00143 CoinCopyN (Ub (), ncols, oub);
00144
00145
00146 CoinCopyN (chg_bds, ncols, f_chg);
00147
00148
00149
00150 break;
00151 }
00152
00153 if ( (direction && ((fb < inner) || (fb > outer))) ||
00154 (!direction && ((fb > inner) || (fb < outer))))
00155 fb = 0.5 * (inner + outer);
00156 }
00157
00158 if (direction) {
00159 Lb (index) = intvar ? ceil (fb - COUENNE_EPS) : fb;
00160 f_chg [index].setLower (t_chg_bounds::CHANGED);
00161 } else {
00162 Ub (index) = intvar ? floor (fb + COUENNE_EPS) : fb;
00163 f_chg [index].setUpper (t_chg_bounds::CHANGED);
00164 }
00165
00166
00167
00168 if (jnlst_ -> ProduceOutput (Ipopt::J_ERROR, J_BOUNDTIGHTENING)) {
00169 char c1 = direction ? '-' : '>', c2 = direction ? '<' : '-';
00170 printf (" # x%d = %g iter %3d: [%+10g -%c %+10g %c- %+10g] /\\/\\ ",index,xcur,iter,olb[index],c1,fb,c2, oub [index]);
00171 printf (" [%10g,%10g]<%g,%g>=> ",Lb (index),Ub (index),CoinMin(inner,outer),CoinMax(inner,outer));
00172 }
00173
00174 bool
00175 feasible = btCore (f_chg),
00176 betterbds = Lb (objind) > getCutOff () + COUENNE_EPS;
00177
00178 jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING,
00179 " [%10g,%10g] lb = %g {fea=%d,btr=%d} ",
00180 Lb (index), Ub (index), Lb (objind),feasible,betterbds);
00181
00182 if (feasible && !betterbds) {
00183
00184
00185 inner = fb;
00186
00187
00188 CoinCopyN (chg_bds, ncols, f_chg);
00189 CoinCopyN (olb, ncols, Lb ());
00190 CoinCopyN (oub, ncols, Ub ());
00191
00192 } else {
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 if (optimum_ &&
00206 ((!direction &&
00207 (optimum_ [index] >= olb [index]) &&
00208 (optimum_ [index] <= Lb (index) - COUENNE_EPS)) ||
00209 (direction &&
00210 (optimum_ [index] <= oub [index]) &&
00211 (optimum_ [index] >= Ub (index) + COUENNE_EPS)))) {
00212
00213 jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING,
00214 "fake tightening CUTS optimum: x%d=%g in [%g,%g] but not in [%g,%g]\n",
00215 index, olb [index], oub [index], Lb (index), Ub (index));
00216 }
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 if (direction) {
00243
00244 oub [index] = Ub (index) = intvar ? floor (fb + COUENNE_EPS) : fb;
00245 chg_bds [index]. setUpper (t_chg_bounds::CHANGED);
00246
00247 } else {
00248
00249 olb [index] = Lb (index) = intvar ? ceil (fb - COUENNE_EPS) : fb;
00250 chg_bds [index]. setLower (t_chg_bounds::CHANGED);
00251 }
00252
00253 outer = fb;
00254
00255 tightened = true;
00256
00257
00258
00259 CoinCopyN (chg_bds, ncols, f_chg);
00260 CoinCopyN (olb, ncols, Lb ());
00261 CoinCopyN (oub, ncols, Ub ());
00262
00263
00264
00265 if (!(btCore (chg_bds))) {
00266
00267 jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING,
00268 "\n pruned by Probing\n");
00269 return -1;
00270
00271 } else {
00272
00273
00274
00275 CoinCopyN (Lb (), ncols, olb);
00276 CoinCopyN (Ub (), ncols, oub);
00277 }
00278
00279 }
00280
00281
00282
00283
00284
00285
00286
00287
00288 CouNumber
00289 lb = Lb (index),
00290 ub = Ub (index);
00291
00292 if ((!direction && ((inner > ub) || (outer < lb))) ||
00293 ( direction && ((inner < lb) || (outer > ub)))) {
00294
00295
00296
00297 inner = direction ? lb : ub;
00298 outer = direction ? ub : lb;
00299 }
00300
00301 CouNumber diff = fabs (inner - outer);
00302
00303 if (diff <= COUENNE_EPS) break;
00304
00305 if (diff > 1.) {
00306
00307 CouNumber L = log (diff) / log (10.);
00308
00309 if (direction) fb = inner + exp (log (10.) * L/2);
00310 else fb = inner - exp (log (10.) * L/2);
00311
00312 } else fb = (inner + outer)/2;
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, "\n");
00323 }
00324
00325 Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, "\n");
00326 if (tightened)
00327 Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING,
00328 " [x%2d] pruned %s [%g, %g] -- lb = %g cutoff = %g\n",
00329 index,direction?"right":"left",
00330 olb[index],oub[index], Lb (objind), getCutOff ());
00331
00332 return tightened ? 1 : 0;
00333 }