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