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