00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CouennePrecisions.hpp"
00012 #include "CouenneTypes.hpp"
00013 #include "CouenneObject.hpp"
00014
00015 #include "exprMul.hpp"
00016 #include "funtriplets.hpp"
00017 #include "projections.hpp"
00018
00019
00022 CouNumber exprMul::selectBranch (const CouenneObject *obj,
00023 const OsiBranchingInformation *info,
00024 expression *&var,
00025 double * &brpts,
00026 double * &brDist,
00027
00028 int &way) {
00029
00030 if (brDist) {free (brDist); brDist = NULL;}
00031
00032 int xi = arglist_ [0] -> Index (),
00033 yi = arglist_ [1] -> Index (),
00034 wi = obj -> Reference () -> Index ();
00035
00036 assert ((xi >= 0) && (yi >= 0) && (wi >= 0));
00037
00038 CouNumber
00039 x0 = info -> solution_ [xi], y0 = info -> solution_ [yi],
00040 xl = info -> lower_ [xi], yl = info -> lower_ [yi],
00041 xu = info -> upper_ [xi], yu = info -> upper_ [yi];
00042
00043
00044 #ifdef DEBUG
00045 printf (" branch MUL: %g [%g,%g] %g [%g,%g]\n",
00046 x0, xl, xu, y0, yl, yu);
00047 #endif
00048
00049 brpts = (double *) realloc (brpts, sizeof (double));
00050
00051
00052
00053 if (fabs (xu-xl) < COUENNE_EPS) {
00054
00055 if (fabs (yu-yl) < COUENNE_EPS) {
00056
00057 var = NULL;
00058 return 0.;
00059
00060 } else {
00061
00062 var = arglist_ [1];
00063 *brpts = 0.5 * (yl+yu);
00064 brDist = (double *) realloc (brDist, 2 * sizeof (double));
00065
00066 brDist [0] = projectSeg (x0, y0, yl, xl*yl, *brpts, *brpts * xl, 0);
00067 brDist [1] = projectSeg (x0, y0, *brpts, *brpts * xl, yu, xl*yu, 0);
00068
00069
00070 return CoinMin (brDist [0], brDist [1]);
00071 }
00072
00073 } else if (fabs (yu-yl) < COUENNE_EPS) {
00074
00075 var = arglist_ [0];
00076 *brpts = 0.5 * (xl+xu);
00077 brDist = (double *) realloc (brDist, 2 * sizeof (double));
00078
00079 brDist [0] = projectSeg (x0, y0, xl, xl*yl, *brpts, *brpts * yl, 0);
00080 brDist [1] = projectSeg (x0, y0, *brpts, *brpts * yl, xu, xu*yl, 0);
00081
00082
00083 return CoinMin (brDist [0], brDist [1]);
00084 }
00085
00086
00087
00088 if (((var = arglist_ [0]) -> Index() >= 0) && (xl < -COUENNE_INFINITY) && (xu > COUENNE_INFINITY) ||
00089 ((var = arglist_ [1]) -> Index() >= 0) && (yl < -COUENNE_INFINITY) && (yu > COUENNE_INFINITY)) {
00090
00091 *brpts = 0.;
00092 brDist = computeMulBrDist (info, xi, yi, wi, var -> Index (), brpts);
00093 way = (info -> solution_ [var -> Index ()] > *brpts) ? TWO_RIGHT : TWO_LEFT;
00094
00095 return CoinMin (brDist [0], brDist [1]);
00096 }
00097
00098
00099
00100
00101
00102 int ind = -1;
00103
00104 if (xl < -COUENNE_INFINITY)
00105 {ind = xi; *brpts = obj -> midInterval (((x0 < 0.) ? 2 : 0.5) * x0, xl, xu); way = TWO_RIGHT;}
00106
00107 else if (xu > COUENNE_INFINITY)
00108 {ind = xi; *brpts = obj -> midInterval (((x0 > 0.) ? 2 : 0.5) * x0, xl, xu); way = TWO_LEFT;}
00109
00110 else if (yl < -COUENNE_INFINITY)
00111 {ind = yi; *brpts = obj -> midInterval (((y0 < 0.) ? 2 : 0.5) * y0, yl, yu); way = TWO_RIGHT;}
00112
00113 else if (yu > COUENNE_INFINITY)
00114 {ind = yi; *brpts = obj -> midInterval (((y0 > 0.) ? 2 : 0.5) * y0, yl, yu) ;way = TWO_LEFT;}
00115
00116 else {
00117
00118
00119
00120 CouNumber delta = (yu-yl) - (xu-xl);
00121
00122 if (delta > +COUENNE_EPS) ind = yi;
00123 else if (delta < -COUENNE_EPS) ind = xi;
00124 else ind = (CoinDrand48 () < 0.5) ? xi : yi;
00125
00126 CouNumber
00127 pt = info -> solution_ [ind],
00128 lb = info -> lower_ [ind],
00129 ub = info -> upper_ [ind];
00130
00131 if ((lb < -COUENNE_EPS) && (ub > COUENNE_EPS) &&
00132 (-lb/ub >= THRES_ZERO_SYMM) &&
00133 (-ub/lb >= THRES_ZERO_SYMM))
00134
00135 *brpts = 0.;
00136
00137 else switch (obj -> Strategy ()) {
00138 case CouenneObject::MID_INTERVAL: *brpts = obj -> midInterval (pt, lb, ub); break;
00139 case CouenneObject::BALANCED: *brpts = balancedMul (info, (ind == xi) ? 0 : 1, wi); break;
00140 case CouenneObject::MIN_AREA:
00141 default: *brpts = (0.5 * (lb+ub)); break;
00142 }
00143
00144 way = (pt > *brpts) ? TWO_RIGHT : TWO_LEFT;
00145 }
00146
00147 assert (ind >= 0);
00148
00149 var = arglist_ [(ind == xi) ? 0 : 1];
00150
00151 brDist = computeMulBrDist (info, xi, yi, wi, ind, brpts);
00152
00153 #ifdef DEBUG
00154 printf (" MUL: br on x_%d %g [%g,%g] [%g,%g] (%g,%g)\n",
00155 ind, *brpts, xl, xu, yl, yu, x0, y0);
00156 #endif
00157
00158 return CoinMin (brDist [0], brDist [1]);
00159
00160 }
00161
00162
00163
00164 CouNumber exprMul::balancedMul (const OsiBranchingInformation *info, int index, int wind) {
00165
00166
00167
00168 int other;
00169
00170 if (index==0) {
00171 index = arglist_ [0] -> Index ();
00172 other = arglist_ [1] -> Index ();
00173 } else {
00174 index = arglist_ [1] -> Index ();
00175 other = arglist_ [0] -> Index ();
00176 }
00177
00178 assert ((index >= 0) && (other >= 0));
00179
00180 CouNumber
00181 xl = info -> lower_ [index], yl = info -> lower_ [other],
00182 xu = info -> upper_ [index], yu = info -> upper_ [other],
00183 x0 = info -> solution_ [index], y0 = info -> solution_ [other],
00184 w0 = info -> solution_ [wind];
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
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
00243
00244
00245
00246
00247
00248
00249
00250 powertriplet ft (2);
00251
00252
00253
00254
00255 bool above = (w0 > x0*y0);
00256
00257 CouNumber
00258 dx = xu-xl,
00259 dy = yu-yl,
00260 area = dx*dy,
00261 bA = yl*dx - xl*dy,
00262 bB = -yl*dx - xu*dy,
00263 m = 1. / sqrt (area),
00264 qA = -bA / (2*area),
00265 qB = bB / (2*area),
00266 z_opt = above ?
00267 minMaxDelta (&ft, -qA/m, (1-qA)/m):
00268 minMaxDelta (&ft, -qB/m, (1-qB)/m),
00269 t_optA = m*z_opt + qA,
00270 t_optB = m*z_opt + qB;
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 return (w0 > x0*y0) ?
00281 (xl + dx * t_optA) :
00282 (xu - dx * t_optB);
00283 }