00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CouenneProblem.hpp"
00012 #include "CouenneExprConst.hpp"
00013 #include "CouenneExprQuad.hpp"
00014 #include "CouenneExprPow.hpp"
00015 #include "CouenneExprMul.hpp"
00016 #include "CouenneDepGraph.hpp"
00017 #include "CouenneLQelems.hpp"
00018
00019 #include "CoinHelperFunctions.hpp"
00020 #include "CoinFinite.hpp"
00021
00022 using namespace Couenne;
00023
00024 namespace Couenne {
00025 class Domain;
00026 }
00027
00028
00029
00030 struct cmpVar {
00031 bool operator() (const exprVar* v1, const exprVar* v2) const
00032 {return (v1 -> Index () < v2 -> Index ());}
00033 };
00034
00036 exprQuad::exprQuad (CouNumber c0,
00037 std::vector <std::pair <exprVar *, CouNumber> > &lcoeff,
00038 std::vector <quadElem> &qcoeff,
00039 expression **al,
00040 int n):
00041
00042 exprGroup (c0, lcoeff, al, n) {
00043
00044 nqterms_ = 0;
00045
00046 typedef std::map <exprVar *, CouNumber, cmpVar> rowMap;
00047 typedef std::map <exprVar *, rowMap, cmpVar> matrixMap;
00048
00049 matrixMap qMap;
00050
00051 for (std::vector <quadElem>::iterator qel = qcoeff.begin (); qel != qcoeff.end (); ++qel) {
00052
00053 CouNumber coe = qel -> coeff ();
00054
00055 exprVar
00056 *varI = qel -> varI (),
00057 *varJ = qel -> varJ ();
00058
00059 if (varI -> Index () != varJ -> Index ()) {
00060
00061
00062
00063
00064 if (varI -> Index () > varJ -> Index ()) {
00065
00066 exprVar *swap = varJ;
00067 varJ = varI;
00068 varI = swap;
00069 }
00070 }
00071
00072 matrixMap::iterator rowp = qMap.find (varI);
00073
00074 if (rowp == qMap.end ()) {
00075
00076 std::pair <exprVar *, CouNumber> newcell (varJ, coe);
00077 rowMap rmap;
00078 rmap.insert (newcell);
00079
00080 std::pair <exprVar *, rowMap> newrow (varI, rmap);
00081 qMap.insert (newrow);
00082
00083 } else {
00084
00085 rowMap::iterator cell = rowp -> second.find (varJ);
00086
00087 if (cell == rowp -> second.end ()) {
00088
00089 std::pair <exprVar *, CouNumber> newcell (varJ, coe);
00090 rowp -> second.insert (newcell);
00091
00092 } else {
00093
00094 if (fabs (cell -> second += coe) < COUENNE_EPS)
00095
00096 rowp -> second.erase (cell);
00097 }
00098 }
00099 }
00100
00101
00102
00103 for (matrixMap::iterator row = qMap.begin (); row != qMap.end (); ++row) {
00104
00105 sparseQcol line;
00106
00107
00108 if (bounds_.find (row -> first) == bounds_.end ()) {
00109
00110 std::pair <CouNumber, CouNumber> newbound (-COIN_DBL_MAX, COIN_DBL_MAX);
00111 std::pair <exprVar *, std::pair <CouNumber, CouNumber> > newvar (row -> first, newbound);
00112 bounds_.insert (newvar);
00113 }
00114
00115 for (rowMap::iterator cell = row -> second.begin (); cell != row -> second.end (); ++cell) {
00116
00117 line.push_back (std::pair <exprVar *, CouNumber> (*cell));
00118
00119
00120 if (bounds_.find (cell -> first) == bounds_.end ()) {
00121
00122 std::pair <CouNumber, CouNumber> newbound (-COIN_DBL_MAX, COIN_DBL_MAX);
00123 std::pair <exprVar *, std::pair <CouNumber, CouNumber> > newvar (cell -> first, newbound);
00124 bounds_.insert (newvar);
00125 }
00126 }
00127
00128 matrix_.push_back (std::pair <exprVar *, sparseQcol> (row -> first, line));
00129 nqterms_ += (int) (line.size ());
00130 }
00131 }
00132
00133
00135 exprQuad::exprQuad (const exprQuad &src, Domain *d):
00136 exprGroup (src, d),
00137 bounds_ (src.bounds_),
00138 nqterms_ (src.nqterms_) {
00139
00140 for (sparseQ::iterator row = src.matrix_.begin (); row != src.matrix_ . end (); ++row) {
00141
00142 sparseQcol column;
00143
00144 for (sparseQcol::iterator i = row -> second. begin (); i != row -> second. end (); ++i)
00145 column.push_back (std::pair <exprVar *, CouNumber>
00146
00147 (new exprVar (i -> first -> Index (), d), i -> second));
00148
00149 matrix_.push_back (std::pair <exprVar *, sparseQcol>
00150
00151 (new exprVar (row -> first -> Index (), d), column));
00152 }
00153
00155
00156 std::vector
00157 <std::pair <CouNumber, std::vector
00158 <std::pair <exprVar *, CouNumber> > > >::iterator row;
00159
00160 for (row = src.eigen_ . begin ();
00161 row != src.eigen_ . end (); ++row) {
00162
00163 std::vector <std::pair <exprVar *, CouNumber> > eigVec;
00164
00165 for (std::vector <std::pair <exprVar *, CouNumber> >::iterator
00166 i = row -> second. begin ();
00167 i != row -> second. end (); ++i)
00168 eigVec.push_back (std::pair <exprVar *, CouNumber>
00169 (dynamic_cast <exprVar *> (i -> first -> clone (d)), i -> second));
00170
00171 eigen_.push_back (std::pair <CouNumber, std::vector
00172 <std::pair <exprVar *, CouNumber> > > (row -> first, eigVec));
00173 }
00174 }
00175
00176
00178 void exprQuad::print (std::ostream &out, bool descend) const {
00179
00180
00181 if (matrix_.size () > 0)
00182 out << '(';
00183
00184
00185 exprGroup::print (out, descend);
00186
00187 int noperands = 0;
00188
00189 for (size_t n = matrix_.size (), i=0; n--; i++) {
00190
00191
00192 int xind = matrix_ [i].first -> Index ();
00193 const sparseQcol row = matrix_ [i].second;
00194
00195 for (int m = row.size (), j=0; m--; j++) {
00196
00197
00198 if (fabs (row [j]. second - 1.) > COUENNE_EPS) {
00199 if (fabs (row [j]. second + 1.) < COUENNE_EPS) out << "- ";
00200 else {
00201 if (row [j]. second > 0.) out << '+';
00202 out << row [j]. second << "*";
00203 }
00204 } else out << '+';
00205
00206 if (row [j].first -> Index () == xind) {
00207 matrix_ [i]. first -> print (out, descend);
00208 out << "^2";
00209 } else {
00210 matrix_ [i]. first -> print (out, descend);
00211 out << '*';
00212 row [j]. first -> print (out, descend);
00213 }
00214
00215 if (!((noperands + 1) % MAX_ARG_LINE))
00216 out << std::endl;
00217 }
00218 }
00219
00220
00221 if (matrix_.size () > 0)
00222 out << ')';
00223 }
00224
00225
00227 expression *exprQuad::differentiate (int index) {
00228
00229 std::map <exprVar *, CouNumber> lmap;
00230
00231 CouNumber c0 = 0;
00232
00233
00234 for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00235 c0 += el -> second;
00236
00237
00238 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00239
00240 int xind = row -> first -> Index ();
00241
00242 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00243
00244 int yind = col -> first -> Index ();
00245
00246 CouNumber coe = col -> second;
00247 exprVar *var = col -> first;
00248
00249 if (xind == index)
00250 if (yind == index) {var = col -> first; coe *= 2;}
00251 else var = col -> first;
00252 else if (yind == index) var = row -> first;
00253 else continue;
00254
00255 std::map <exprVar *, CouNumber>::iterator i = lmap.find (var);
00256
00257 if (i != lmap.end()) {
00258 if (fabs (i -> second += coe) < COUENNE_EPS)
00259 lmap.erase (i);
00260 } else {
00261 std::pair <exprVar *, CouNumber> npair (var, coe);
00262 lmap.insert (npair);
00263 }
00264 }
00265 }
00266
00267
00268 expression **arglist = new expression * [nargs_ + 1];
00269 int nargs = 0;
00270
00271 for (int i = 0; i < nargs_; i++)
00272 if (arglist_ [i] -> dependsOn (index))
00273 arglist [nargs++] = arglist_ [i] -> differentiate (index);
00274
00275
00276
00277
00278 if (lmap.empty ()) {
00279
00280
00281 if (!nargs) {
00282 delete arglist;
00283 return new exprConst (c0);
00284 }
00285
00286 if (fabs (c0) > COUENNE_EPS)
00287 arglist [nargs++] = new exprConst (c0);
00288
00289 return new exprSum (arglist, nargs);
00290 }
00291
00292 lincoeff coe;
00293
00294 for (std::map <exprVar *, CouNumber>::iterator i = lmap.begin (); i != lmap.end (); ++i)
00295 coe.push_back (std::pair <exprVar *, CouNumber> (i -> first, i -> second));
00296
00297 return new exprGroup (c0, coe, arglist, nargs);
00298 }
00299
00300
00302
00303 int exprQuad::compare (exprQuad &e) {
00304
00305 int sum = exprGroup::compare (e);
00306
00307 if (sum != 0)
00308 return sum;
00309
00310 if (matrix_.size() < e.matrix_.size()) return -1;
00311 if (matrix_.size() > e.matrix_.size()) return 1;
00312
00313 for (sparseQ::iterator
00314 row1 = matrix_.begin (),
00315 row2 = e.matrix_.begin ();
00316 row1 != matrix_.end ();
00317 ++row1, ++row2) {
00318
00319 if (row1 -> first -> Index () < row2 -> first -> Index ()) return -1;
00320 if (row1 -> first -> Index () > row2 -> first -> Index ()) return 1;
00321
00322 if (row1 -> second.size () < row2 -> second.size ()) return -1;
00323 if (row1 -> second.size () > row2 -> second.size ()) return 1;
00324
00325
00326
00327
00328
00329 for (sparseQcol::iterator
00330 col1 = row1 -> second.begin (),
00331 col2 = row2 -> second.begin ();
00332 col1 != row1 -> second.end ();
00333 ++col1, ++col2) {
00334
00335 if (col1 -> first -> Index () < col2 -> first -> Index ()) return -1;
00336 if (col1 -> first -> Index () > col2 -> first -> Index ()) return 1;
00337
00338 if (col1 -> second < col2 -> second - COUENNE_EPS) return -1;
00339 if (col1 -> second > col2 -> second + COUENNE_EPS) return 1;
00340 }
00341 }
00342
00343 return 0;
00344 }
00345
00346
00348
00349 int exprQuad::rank () {
00350
00351 int maxrank = exprGroup::rank ();
00352
00353 if (maxrank < 0)
00354 maxrank = 0;
00355
00356 int r;
00357
00358 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00359
00360 if ((r = row -> first -> rank ()) > maxrank) maxrank = r;
00361
00362 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
00363 if ((r = col -> first -> rank ()) > maxrank) maxrank = r;
00364 }
00365
00366 return maxrank;
00367 }
00368
00369
00371 void exprQuad::fillDepSet (std::set <DepNode *, compNode> *dep, DepGraph *g) {
00372
00373 exprGroup::fillDepSet (dep, g);
00374
00375 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00376
00377 dep -> insert (g -> lookup (row -> first -> Index ()));
00378
00379 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
00380 dep -> insert (g -> lookup (col -> first -> Index ()));
00381 }
00382 }
00383
00384
00387 int exprQuad::DepList (std::set <int> &deplist,
00388 enum dig_type type) {
00389
00390 int deps = exprGroup::DepList (deplist, type);
00391
00392 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00393 deps += row -> first -> DepList (deplist, type);
00394 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
00395 deps += col -> first -> DepList (deplist, type);
00396 }
00397
00398 return deps;
00399 }
00400
00401
00403 bool exprQuad::isInteger () {
00404
00405 if (!(exprGroup::isInteger ()))
00406 return false;
00407
00408 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00409
00410 bool intI = row -> first -> isInteger ();
00411
00412 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00413
00414 CouNumber coe = col -> second;
00415
00416 bool
00417 intCoe = ::isInteger (coe),
00418 intJ = row -> first -> isInteger ();
00419
00420 if (intI && intJ && intCoe)
00421 continue;
00422
00423 if (!intCoe
00424 && row -> first -> isFixed ()
00425 && col -> first -> isFixed ()
00426 && ::isInteger (coe *
00427 row -> first -> lb () *
00428 col -> first -> lb ()))
00429 continue;
00430
00431 if (!intI && (row -> first -> isFixed ()) && ::isInteger ((*(row -> first)) ())) continue;
00432 if (!intJ && (col -> first -> isFixed ()) && ::isInteger ((*(col -> first)) ())) continue;
00433
00434
00435
00436
00437
00438 return false;
00439 }
00440 }
00441
00442 return true;
00443 }
00444
00445
00447 void exprQuad::replace (exprVar *x, exprVar *w) {
00448
00449 exprGroup::replace (x, w);
00450 int xind = x -> Index ();
00451 int wind = w -> Index ();
00452
00453 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00454
00455 exprVar * &vr = row -> first;
00456 if ((vr -> Index () == xind)) {
00457
00458
00459
00460 vr = w;
00461 }
00462
00463 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00464
00465 exprVar * &vc = col -> first;
00466 if ((vc -> Index () == wind)) {
00467
00468
00469
00470 vc = w;
00471 }
00472 }
00473 }
00474 }
00475
00476
00478 void exprQuad::realign (const CouenneProblem *p) {
00479
00480 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00481
00482 exprVar * &vr = row -> first;
00483 int indVar;
00484
00485
00486
00487 if (((vr -> Type () == VAR) ||
00488 (vr -> Type () == AUX)) &&
00489 (vr -> Original () != p -> Var (indVar = vr -> Index ()))) {
00490
00491 expression *trash = vr;
00492 row -> first = p -> Var (indVar);
00493 delete trash;
00494 }
00495
00496
00497
00498 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00499
00500 exprVar * &vc = col -> first;
00501 int indVar;
00502
00503
00504
00505 if (((vc -> Type () == VAR) ||
00506 (vc -> Type () == AUX)) &&
00507 (vc -> Original () != p -> Var (indVar = vc -> Index ()))) {
00508
00509 expression *trash = vc;
00510 col -> first = p -> Var (indVar);
00511 delete trash;
00512 }
00513 }
00514 }
00515 }
00516
00517
00519 void exprQuad::closestFeasible (expression *varind,
00520 expression *vardep,
00521 CouNumber &left,
00522 CouNumber &right) const {
00523
00524 fprintf (stderr, "exprQuad::closestFeasible() not available for VT\n");
00525 exit (-1);
00526 }
00527
00528
00530 CouNumber exprQuad::gradientNorm (const double *x) {
00531
00532 CouNumber grad = 0.;
00533
00534 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00535
00536 CouNumber gradEl = 0.;
00537 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
00538 gradEl += col -> second * x [col -> first -> Index ()];
00539
00540 grad += gradEl * gradEl;
00541 }
00542
00543 return sqrt (grad);
00544 }
00545
00547 expression *exprQuad::simplify () {
00548 exprOp::simplify ();
00549 return NULL;
00550 }
00551