00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CouenneProblem.hpp"
00012
00013 #include "CoinHelperFunctions.hpp"
00014
00015 #include "exprClone.hpp"
00016 #include "exprSum.hpp"
00017 #include "exprMul.hpp"
00018 #include "exprGroup.hpp"
00019 #include "exprQuad.hpp"
00020 #include "lqelems.hpp"
00021
00022
00023
00027 void elementBreak (expression *, int &, CouNumber &);
00028
00029
00032
00033 int CouenneProblem::splitAux (CouNumber rhs, expression *body, expression *&rest, bool *wentAux) {
00034
00035 int auxInd = -1,
00036 code = body -> code ();
00037
00038 expression **alist = body -> ArgList ();
00039
00040 #ifdef DEBUG
00041 printf ("|||||||||| Splitting "); body -> print (); printf ("\n");
00042 #endif
00043
00044 switch (code) {
00045
00046
00048
00049 case COU_EXPRSUB: {
00050
00051
00052
00053 int pos = 0;
00054 CouNumber coeff = 1;
00055
00056 auxInd = (*alist) -> Index ();
00057
00058 if (auxInd < 0)
00059 elementBreak (*alist, auxInd, coeff);
00060
00061 if ((auxInd < 0) ||
00062 wentAux [auxInd] ||
00063 (alist [1] -> dependsOn (auxInd, TAG_AND_RECURSIVE) >= 1)) {
00064
00065 auxInd = (alist [1]) -> Index ();
00066
00067 if (auxInd < 0)
00068 elementBreak (alist [1], auxInd, coeff);
00069 else coeff = 1;
00070
00071 if ((auxInd < 0) ||
00072 wentAux [auxInd] ||
00073 (alist [0] -> dependsOn (auxInd, TAG_AND_RECURSIVE) >= 1))
00074
00075 return -1;
00076
00077 pos = 1;
00078 }
00079
00081
00082
00083 expression *clone = alist [1 - pos] -> clone (&domain_);
00084
00085 expression *auxdef =
00086 (fabs (coeff - 1) < COUENNE_EPS) ?
00087 (clone) :
00088 new exprMul (new exprConst (coeff), clone);
00089
00090 rest = (fabs (rhs) < COUENNE_EPS) ?
00091 (auxdef) :
00092 (new exprSum (auxdef, new exprConst ((pos==1) ? -rhs : rhs)));
00093
00094 } break;
00095
00097
00098 case COU_EXPRQUAD:
00099 case COU_EXPRGROUP:
00100 case COU_EXPRSUM: {
00101
00102
00103
00104
00105
00106
00107
00108
00109 int maxindex = -1, nlin = 0, which = 1;
00110 CouNumber c0 = 0., auxcoe = 1;
00111 bool which_was_set = false;
00112
00113
00114
00115 if (code != COU_EXPRSUM) {
00116
00117 exprGroup *egBody = dynamic_cast <exprGroup *> (body -> isaCopy () ?
00118 body -> Copy () :
00119 body);
00120 exprGroup::lincoeff &lcoe = egBody -> lcoeff ();
00121
00122
00123
00124 c0 = egBody -> getc0 ();
00125
00126 for (int i=0, n = lcoe.size (); n--; i++, nlin++) {
00127
00128 int j = lcoe [i]. first -> Index ();
00129
00130
00131
00132
00133 if ((!(lcoe [i].first -> isInteger ()) || (which==1)) &&
00134 (j > maxindex) &&
00135 !(wentAux [j]) &&
00136 (fabs (lcoe [i]. second) > COUENNE_EPS)) {
00137
00138
00139
00140
00141 exprVar *saveVar = lcoe [i].first;
00142
00143 lcoe [i].first = new exprVar (-1);
00144
00145 if (body -> dependsOn (j, TAG_AND_RECURSIVE) == 0) {
00146
00147 #ifdef DEBUG
00148 printf ("body does not depend on x%d: ", j);
00149 body -> print ();
00150 printf ("\n");
00151 #endif
00152
00153 which = - nlin - 1;
00154 auxcoe = lcoe [i]. second;
00155 maxindex = j;
00156 }
00157
00158 delete lcoe [i].first;
00159 lcoe [i].first = saveVar;
00160 }
00161 }
00162 }
00163
00164 if (which != 1)
00165 which_was_set = true;
00166 else which = -1;
00167
00168
00169
00170 for (int i = body -> nArgs (); i--;) {
00171
00172 CouNumber coeff = 1;
00173 int index = alist [i] -> Index ();
00174
00175 if (index < 0)
00176 elementBreak (alist [i], index, coeff);
00177
00178 if ((index > maxindex) &&
00179 !(wentAux [index]) &&
00180 (fabs (coeff) > COUENNE_EPS)) {
00181
00182
00183
00184 expression *cut = alist [i];
00185 alist [i] = new exprConst (0.);
00186
00187
00188
00189 if (body -> dependsOn (index, TAG_AND_RECURSIVE) == 0) {
00190
00191 maxindex = index;
00192 which = i;
00193 auxcoe = coeff;
00194 }
00195
00196 delete alist [i];
00197 alist [i] = cut;
00198 }
00199 }
00200
00201 if (!which_was_set && (which == -1))
00202 return -1;
00203
00205
00206 if (maxindex < 0) break;
00207
00208
00209
00210
00211
00212 int nargs = body -> nArgs ();
00213 expression **newarglist;
00214
00215 #ifdef DEBUG
00216 printf (" [[ind %d, coe %.1g, wh %d, nargs %d, nlin %d]] ",
00217 maxindex, auxcoe, which, nargs, nlin);
00218 #endif
00219
00220 if (nargs > 0) {
00221
00222 int j, mid = (which < 0) ? nargs : which;
00223
00224 newarglist = new expression * [nargs + 1];
00225
00226 for (j=0; j<mid; j++) newarglist [j] = alist [j] -> clone (&domain_);
00227 for (j++; j<nargs; j++) newarglist [j-1] = alist [j] -> clone (&domain_);
00228
00229
00230
00231
00232 } else {
00233
00234 nargs++;
00235 newarglist = new expression *;
00236 *newarglist = new exprConst (0.);
00237 }
00238
00239
00240
00241 int *linind2 = NULL;
00242 CouNumber *lincoe2 = NULL;
00243
00244
00245 int *qindI = NULL,
00246 *qindJ = NULL;
00247 CouNumber *qcoe = NULL;
00248
00249 if (nlin > 0) {
00250
00251 exprGroup *egBody = dynamic_cast <exprGroup *> (body -> isaCopy () ?
00252 body -> Copy () :
00253 body);
00254 exprGroup::lincoeff &lcoe = egBody -> lcoeff ();
00255
00256 int mid = (which >= 0) ? nlin : - which - 1;
00257
00258 linind2 = new int [nlin + 1];
00259 lincoe2 = new CouNumber [nlin + 1];
00260
00261 register int j;
00262
00263 #ifdef DEBUG
00264
00265
00266 #endif
00267
00268 CouNumber divider = -1. / auxcoe;
00269
00270 for (j=0; j<mid; j++){linind2[j] =lcoe[j].first->Index();lincoe2[j] =divider*lcoe[j].second;}
00271 for (j++; j<nlin; j++){linind2[j-1]=lcoe[j].first->Index();lincoe2[j-1]=divider*lcoe[j].second;}
00272
00273 linind2 [j-1] = -1;
00274
00275 #ifdef DEBUG
00276 for (j=0; j<mid; j++) printf ("<%g x%d> ", lincoe2 [j], linind2 [j]);
00277 for (j++; j<nlin; j++) printf ("<%g x%d> ", lincoe2 [j-1], linind2 [j-1]);
00278 #endif
00279
00280 if (code == COU_EXPRQUAD) {
00281
00282 exprQuad *eq = dynamic_cast <exprQuad *> (body -> isaCopy () ?
00283 body -> Copy () :
00284 body);
00285
00286 int nqt = eq -> getnQTerms (), j=0;
00287
00288 qindI = new int [nqt];
00289 qindJ = new int [nqt];
00290 qcoe = new CouNumber [nqt];
00291
00292 exprQuad::sparseQ &M = eq -> getQ ();
00293
00294 for (exprQuad::sparseQ::iterator row = M.begin ();
00295 row != M.end (); ++row) {
00296
00297 int xind = row -> first -> Index ();
00298
00299 for (exprQuad::sparseQcol::iterator col = row -> second.begin ();
00300 col != row -> second.end (); ++col, ++j) {
00301
00302 qindI [j] = xind;
00303 qindJ [j] = col -> first -> Index ();
00304 qcoe [j] = divider * col -> second;
00305 }
00306 }
00307 }
00308
00309
00310
00311 }
00312
00313
00314 if (which >= 0) --nargs;
00315 else --nlin;
00316
00317 #ifdef DEBUG
00318 printf ("\n::: auxcoe %g, rhs %g, lin %d, nl %d\n", auxcoe, rhs, nlin, nargs);
00319 #endif
00320
00321
00322
00323
00324 if ((code == COU_EXPRQUAD) ||
00325 (code == COU_EXPRGROUP) && (nlin > 0)) {
00326
00327
00328
00329
00330
00331
00332
00333
00334 if (fabs (auxcoe + 1) < COUENNE_EPS) {
00335
00336
00337 exprGroup::lincoeff lcoeff;
00338 indcoe2vector (linind2, lincoe2, lcoeff);
00339
00340 if (code == COU_EXPRGROUP)
00341 rest = new exprGroup (c0-rhs, lcoeff, newarglist, nargs);
00342
00343 else {
00344 std::vector <quadElem> qcoeff;
00345 indcoe2vector (qindI, qindJ, qcoe, qcoeff);
00346 rest = new exprQuad (c0-rhs, lcoeff, qcoeff, newarglist, nargs);
00347 }
00348 }
00349 else {
00350
00351 expression **mullist = new expression * [1];
00352
00353
00354 if ((nargs <= 1) && ((*newarglist) -> Linearity () <= CONSTANT)) {
00355 *mullist = new exprConst ((*newarglist) -> Value ());
00356
00357 delete [] newarglist;
00358 }
00359 else if ((nargs <= 1) &&
00360 ((*newarglist) -> code () == COU_EXPROPP) &&
00361 (fabs (auxcoe - 1) < COUENNE_EPS))
00362
00363 *mullist = (*newarglist) -> Argument ();
00364 else {
00365
00366 if (nargs <= 1)
00367 *mullist = new exprMul (new exprConst (-1. / auxcoe),
00368 *newarglist);
00369 else
00370 *mullist = new exprMul (new exprConst (-1. / auxcoe),
00371 new exprSum (newarglist, nargs));
00372 }
00373
00374 std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
00375 indcoe2vector (linind2, lincoe2, lcoeff);
00376
00377
00378 if (code == COU_EXPRGROUP)
00379 rest = new exprGroup ((rhs - c0) / auxcoe, lcoeff, mullist,1);
00380
00381 else {
00382 std::vector <quadElem> qcoeff;
00383 indcoe2vector (qindI, qindJ, qcoe, qcoeff);
00384 rest = new exprQuad ((rhs - c0) / auxcoe, lcoeff, qcoeff, mullist,1);
00385 }
00386 }
00387 }
00388 else {
00389
00390 if (fabs (rhs) > COUENNE_EPS) {
00391
00392 if ((nargs == 1) && ((*newarglist) -> Type () == CONST)) {
00393
00394 CouNumber val = (*newarglist) -> Value () - rhs;
00395 delete *newarglist;
00396 *newarglist = new exprConst (val);
00397 }
00398 else newarglist [nargs++] = new exprConst (-rhs);
00399 }
00400
00401
00402
00403 expression *auxDef;
00404
00405 if (nargs==1) {
00406 auxDef = *newarglist;
00407 delete [] newarglist;
00408 } else auxDef = new exprSum (newarglist, nargs);
00409
00410 if (fabs (auxcoe + 1) < COUENNE_EPS)
00411 rest = auxDef;
00412 else if ((fabs (auxcoe - 1) < COUENNE_EPS) &&
00413 (auxDef -> code () == COU_EXPROPP))
00414 rest = auxDef -> Argument ();
00415 else
00416
00417 rest = new exprMul (new exprConst (-1./auxcoe),
00418 ((auxDef -> Type () == AUX) ||
00419 (auxDef -> Type () == VAR)) ? new exprClone (auxDef) : auxDef);
00420 }
00421
00422 if (linind2) {
00423 delete [] linind2;
00424 delete [] lincoe2;
00425 }
00426
00427 #ifdef DEBUG
00428 printf ("gotten "); rest -> print (); printf ("\n");
00429 #endif
00430
00431 auxInd = maxindex;
00432
00433 } break;
00434
00435 default: break;
00436 }
00437
00438 if ((auxInd < 0) || (wentAux [auxInd]))
00439 return -1;
00440
00441
00442
00443
00444
00445
00446 #ifdef DEBUG
00447 printf ("standardize rest (2nd level) "); fflush (stdout);
00448 rest -> print (); printf (" [body = ");
00449 body -> print (); printf ("]");
00450 #endif
00451
00452
00453
00454 exprAux *aux = rest -> standardize (this, false);
00455
00456 #ifdef DEBUG
00457 printf (" {body = "); body -> print (); printf ("} ");
00458 if (aux) {
00459 printf (" --- aux = ");
00460 aux -> print ();
00461 printf (" := ");
00462 aux -> Image () -> print ();
00463 }
00464 #endif
00465
00466 if (aux) {
00467
00468
00469 rest = aux -> Image ();
00470 aux -> Image (NULL);
00471 delete aux;
00472 }
00473
00474 #ifdef DEBUG
00475 printf (" ==> "); rest -> print ();
00476 printf (" and "); fflush (stdout);
00477 rest -> print (); printf ("\n");
00478 #endif
00479
00480 return auxInd;
00481 }