1 /* $Id$
2 *
3 * Name: exprQuad.cpp
4 * Author: Pietro Belotti
5 * Purpose: implementation of some methods for exprQuad
6 *
7 * (C) Carnegie-Mellon University, 2006-08.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CouenneProblem.hpp"
12 #include "CouenneExprConst.hpp"
13 #include "CouenneExprQuad.hpp"
14 #include "CouenneExprPow.hpp"
15 #include "CouenneExprMul.hpp"
16 #include "CouenneDepGraph.hpp"
17 #include "CouenneLQelems.hpp"
18
19 #include "CoinHelperFunctions.hpp"
20 #include "CoinFinite.hpp"
21
22 using namespace Couenne;
23
24 namespace Couenne {
25 class Domain;
26 }
27
28 //#define DEBUG
29
30 struct cmpVar {
operator ()cmpVar31 bool operator() (const exprVar* v1, const exprVar* v2) const
32 {return (v1 -> Index () < v2 -> Index ());}
33 };
34
35 /// Constructor
exprQuad(CouNumber c0,std::vector<std::pair<exprVar *,CouNumber>> & lcoeff,std::vector<quadElem> & qcoeff,expression ** al,int n)36 exprQuad::exprQuad (CouNumber c0,
37 std::vector <std::pair <exprVar *, CouNumber> > &lcoeff,
38 std::vector <quadElem> &qcoeff,
39 expression **al,
40 int n):
41
42 exprGroup (c0, lcoeff, al, n) {
43
44 nqterms_ = 0;
45
46 typedef std::map <exprVar *, CouNumber, cmpVar> rowMap;
47 typedef std::map <exprVar *, rowMap, cmpVar> matrixMap;
48
49 matrixMap qMap;
50
51 for (std::vector <quadElem>::iterator qel = qcoeff.begin (); qel != qcoeff.end (); ++qel) {
52
53 CouNumber coe = qel -> coeff ();
54
55 exprVar
56 *varI = qel -> varI (),
57 *varJ = qel -> varJ ();
58
59 if (varI -> Index () != varJ -> Index ()) {
60
61 //coe /= 2.;
62
63 // pick smaller index as row reference
64 if (varI -> Index () > varJ -> Index ()) {
65
66 exprVar *swap = varJ;
67 varJ = varI;
68 varI = swap;
69 }
70 }
71
72 matrixMap::iterator rowp = qMap.find (varI);
73
74 if (rowp == qMap.end ()) { // add new row
75
76 std::pair <exprVar *, CouNumber> newcell (varJ, coe);
77 rowMap rmap;
78 rmap.insert (newcell);
79
80 std::pair <exprVar *, rowMap> newrow (varI, rmap);
81 qMap.insert (newrow);
82
83 } else { // insert element into row
84
85 rowMap::iterator cell = rowp -> second.find (varJ);
86
87 if (cell == rowp -> second.end ()) { // normal case, add entry
88
89 std::pair <exprVar *, CouNumber> newcell (varJ, coe);
90 rowp -> second.insert (newcell);
91
92 } else { // strange, but add coefficient
93
94 if (fabs (cell -> second += coe) < COUENNE_EPS)
95 // eliminate element of map if null coefficient
96 rowp -> second.erase (cell);
97 }
98 }
99 }
100
101 // transform maps into vectors
102
103 for (matrixMap::iterator row = qMap.begin (); row != qMap.end (); ++row) {
104
105 sparseQcol line;
106
107 // insert first element in bound map
108 if (bounds_.find (row -> first) == bounds_.end ()) {
109
110 std::pair <CouNumber, CouNumber> newbound (-COIN_DBL_MAX, COIN_DBL_MAX);
111 std::pair <exprVar *, std::pair <CouNumber, CouNumber> > newvar (row -> first, newbound);
112 bounds_.insert (newvar);
113 }
114
115 for (rowMap::iterator cell = row -> second.begin (); cell != row -> second.end (); ++cell) {
116
117 line.push_back (std::pair <exprVar *, CouNumber> (*cell));
118
119 // insert second element in bound map
120 if (bounds_.find (cell -> first) == bounds_.end ()) {
121
122 std::pair <CouNumber, CouNumber> newbound (-COIN_DBL_MAX, COIN_DBL_MAX);
123 std::pair <exprVar *, std::pair <CouNumber, CouNumber> > newvar (cell -> first, newbound);
124 bounds_.insert (newvar);
125 }
126 }
127
128 matrix_.push_back (std::pair <exprVar *, sparseQcol> (row -> first, line));
129 nqterms_ += (int) (line.size ());
130 }
131 }
132
133
134 /// copy constructor
exprQuad(const exprQuad & src,Domain * d)135 exprQuad::exprQuad (const exprQuad &src, Domain *d):
136 exprGroup (src, d),
137 bounds_ (src.bounds_),
138 nqterms_ (src.nqterms_) {
139
140 for (sparseQ::iterator row = src.matrix_.begin (); row != src.matrix_ . end (); ++row) {
141
142 sparseQcol column;
143
144 for (sparseQcol::iterator i = row -> second. begin (); i != row -> second. end (); ++i)
145 column.push_back (std::pair <exprVar *, CouNumber>
146 //(dynamic_cast <exprVar *> (i -> first -> clone (d)), i -> second));
147 (new exprVar (i -> first -> Index (), d), i -> second));
148
149 matrix_.push_back (std::pair <exprVar *, sparseQcol>
150 //dynamic_cast <exprVar *> (row -> first -> clone (d)), column));
151 (new exprVar (row -> first -> Index (), d), column));
152 }
153
154 //////////////////////////////////////////////////////////////////////////////
155
156 std::vector
157 <std::pair <CouNumber, std::vector
158 <std::pair <exprVar *, CouNumber> > > >::iterator row;
159
160 for (row = src.eigen_ . begin ();
161 row != src.eigen_ . end (); ++row) {
162
163 std::vector <std::pair <exprVar *, CouNumber> > eigVec;
164
165 for (std::vector <std::pair <exprVar *, CouNumber> >::iterator
166 i = row -> second. begin ();
167 i != row -> second. end (); ++i)
168 eigVec.push_back (std::pair <exprVar *, CouNumber>
169 (dynamic_cast <exprVar *> (i -> first -> clone (d)), i -> second));
170
171 eigen_.push_back (std::pair <CouNumber, std::vector
172 <std::pair <exprVar *, CouNumber> > > (row -> first, eigVec));
173 }
174 }
175
176
177 /// I/O
print(std::ostream & out,bool descend) const178 void exprQuad::print (std::ostream &out, bool descend) const {
179
180 //if (code () == COU_EXPRQUAD)
181 if (matrix_.size () > 0)
182 out << '(';
183
184 // print linear and nonquadratic part
185 exprGroup::print (out, descend);
186
187 int noperands = 0;
188
189 for (size_t n = matrix_.size (), i=0; n--; i++) {
190 //sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
191
192 int xind = matrix_ [i].first -> Index ();
193 const sparseQcol row = matrix_ [i].second;
194
195 for (int m = row.size (), j=0; m--; j++) {
196 //sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
197
198 if (fabs (row [j]. second - 1.) > COUENNE_EPS) {
199 if (fabs (row [j]. second + 1.) < COUENNE_EPS) out << "- ";
200 else {
201 if (row [j]. second > 0.) out << '+';
202 out << row [j]. second << "*";
203 }
204 } else out << '+';
205
206 if (row [j].first -> Index () == xind) {
207 matrix_ [i]. first -> print (out, descend);
208 out << "^2";
209 } else {
210 matrix_ [i]. first -> print (out, descend);
211 out << '*';
212 row [j]. first -> print (out, descend);
213 }
214
215 if (!((noperands + 1) % MAX_ARG_LINE))
216 out << std::endl;
217 }
218 }
219
220 //if (code () == COU_EXPRGROUP)
221 if (matrix_.size () > 0)
222 out << ')';
223 }
224
225
226 /// differentiation
differentiate(int index)227 expression *exprQuad::differentiate (int index) {
228
229 std::map <exprVar *, CouNumber> lmap;
230
231 CouNumber c0 = 0;
232
233 // derive linear part (obtain constant)
234 for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
235 c0 += el -> second;
236
237 // derive quadratic part (obtain linear part)
238 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
239
240 int xind = row -> first -> Index ();
241
242 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
243
244 int yind = col -> first -> Index ();
245
246 CouNumber coe = col -> second;
247 exprVar *var = col -> first;
248
249 if (xind == index)
250 if (yind == index) {var = col -> first; coe *= 2;}
251 else var = col -> first;
252 else if (yind == index) var = row -> first;
253 else continue;
254
255 std::map <exprVar *, CouNumber>::iterator i = lmap.find (var);
256
257 if (i != lmap.end()) {
258 if (fabs (i -> second += coe) < COUENNE_EPS)
259 lmap.erase (i);
260 } else {
261 std::pair <exprVar *, CouNumber> npair (var, coe);
262 lmap.insert (npair);
263 }
264 }
265 }
266
267 // derive nonlinear sum
268 expression **arglist = new expression * [nargs_ + 1];
269 int nargs = 0;
270
271 for (int i = 0; i < nargs_; i++)
272 if (arglist_ [i] -> dependsOn (index))
273 arglist [nargs++] = arglist_ [i] -> differentiate (index);
274
275 // special cases
276
277 // 1) no linear part
278 if (lmap.empty ()) {
279
280 // and no nonlinear part either
281 if (!nargs) {
282 delete arglist;
283 return new exprConst (c0);
284 }
285
286 if (fabs (c0) > COUENNE_EPS)
287 arglist [nargs++] = new exprConst (c0);
288
289 return new exprSum (arglist, nargs);
290 }
291
292 lincoeff coe;
293
294 for (std::map <exprVar *, CouNumber>::iterator i = lmap.begin (); i != lmap.end (); ++i)
295 coe.push_back (std::pair <exprVar *, CouNumber> (i -> first, i -> second));
296
297 return new exprGroup (c0, coe, arglist, nargs);
298 }
299
300
301 /// compare quadratic terms
302
compare(exprQuad & e)303 int exprQuad::compare (exprQuad &e) {
304
305 int sum = exprGroup::compare (e);
306
307 if (sum != 0)
308 return sum;
309
310 if (matrix_.size() < e.matrix_.size()) return -1;
311 if (matrix_.size() > e.matrix_.size()) return 1;
312
313 for (sparseQ::iterator
314 row1 = matrix_.begin (),
315 row2 = e.matrix_.begin ();
316 row1 != matrix_.end ();
317 ++row1, ++row2) {
318
319 if (row1 -> first -> Index () < row2 -> first -> Index ()) return -1;
320 if (row1 -> first -> Index () > row2 -> first -> Index ()) return 1;
321
322 if (row1 -> second.size () < row2 -> second.size ()) return -1;
323 if (row1 -> second.size () > row2 -> second.size ()) return 1;
324
325 // if (matrix_.size() > e.matrix_.size()) return 1;
326 // int xind = row -> first -> Index ();
327 // CouNumber x = (*(row -> first)) ();
328
329 for (sparseQcol::iterator
330 col1 = row1 -> second.begin (),
331 col2 = row2 -> second.begin ();
332 col1 != row1 -> second.end ();
333 ++col1, ++col2) {
334
335 if (col1 -> first -> Index () < col2 -> first -> Index ()) return -1;
336 if (col1 -> first -> Index () > col2 -> first -> Index ()) return 1;
337
338 if (col1 -> second < col2 -> second - COUENNE_EPS) return -1;
339 if (col1 -> second > col2 -> second + COUENNE_EPS) return 1;
340 }
341 }
342
343 return 0;
344 }
345
346
347 /// used in rank-based branching variable choice
348
rank()349 int exprQuad::rank () {
350
351 int maxrank = exprGroup::rank ();
352
353 if (maxrank < 0)
354 maxrank = 0;
355
356 int r;
357
358 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
359
360 if ((r = row -> first -> rank ()) > maxrank) maxrank = r;
361
362 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
363 if ((r = col -> first -> rank ()) > maxrank) maxrank = r;
364 }
365
366 return maxrank;
367 }
368
369
370 /// update dependence set with index of this variable
fillDepSet(std::set<DepNode *,compNode> * dep,DepGraph * g)371 void exprQuad::fillDepSet (std::set <DepNode *, compNode> *dep, DepGraph *g) {
372
373 exprGroup::fillDepSet (dep, g);
374
375 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
376
377 dep -> insert (g -> lookup (row -> first -> Index ()));
378
379 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
380 dep -> insert (g -> lookup (col -> first -> Index ()));
381 }
382 }
383
384
385 /// fill in the set with all indices of variables appearing in the
386 /// expression
DepList(std::set<int> & deplist,enum dig_type type)387 int exprQuad::DepList (std::set <int> &deplist,
388 enum dig_type type) {
389
390 int deps = exprGroup::DepList (deplist, type);
391
392 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
393 deps += row -> first -> DepList (deplist, type);
394 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
395 deps += col -> first -> DepList (deplist, type);
396 }
397
398 return deps;
399 }
400
401
402 /// is this quadratic expression integer?
isInteger()403 bool exprQuad::isInteger () {
404
405 if (!(exprGroup::isInteger ()))
406 return false;
407
408 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
409
410 bool intI = row -> first -> isInteger ();
411
412 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
413
414 CouNumber coe = col -> second;
415
416 bool
417 intCoe = ::isInteger (coe),
418 intJ = row -> first -> isInteger ();
419
420 if (intI && intJ && intCoe)
421 continue;
422
423 if (!intCoe // coefficient fractional, check all is fixed and product is integer
424 && row -> first -> isFixed ()
425 && col -> first -> isFixed ()
426 && ::isInteger (coe *
427 row -> first -> lb () *
428 col -> first -> lb ()))
429 continue;
430
431 if (!intI && (row -> first -> isFixed ()) && ::isInteger ((*(row -> first)) ())) continue;
432 if (!intJ && (col -> first -> isFixed ()) && ::isInteger ((*(col -> first)) ())) continue;
433
434 //if (!intI && !intJ && intCoe) ; // check x y fixed int
435 //if (!intI && intJ && intCoe) ; // check x fixed int
436 //if ( intI && !intJ && intCoe) ; // check y fixed int
437
438 return false;
439 }
440 }
441
442 return true;
443 }
444
445
446 /// replace variable x with new (aux) w
replace(exprVar * x,exprVar * w)447 void exprQuad::replace (exprVar *x, exprVar *w) {
448
449 exprGroup::replace (x, w);
450 int xind = x -> Index ();
451 int wind = w -> Index ();
452
453 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
454
455 exprVar * &vr = row -> first;
456 if ((vr -> Index () == xind)) {
457
458 //fprintf (stderr, "Didn't fix exprQuad::replace() yet");
459 //exit (-1);
460 vr = w;
461 }
462
463 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
464
465 exprVar * &vc = col -> first;
466 if ((vc -> Index () == wind)) {
467
468 //fprintf (stderr, "Didn't fix exprQuad::replace() yet");
469 //exit (-1);
470 vc = w;
471 }
472 }
473 }
474 }
475
476
477 /// Redirect variables to proper variable vector
realign(const CouenneProblem * p)478 void exprQuad::realign (const CouenneProblem *p) {
479
480 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
481
482 exprVar * &vr = row -> first;
483 int indVar;
484
485 // substitute variable representing this row with its newest version
486
487 if (((vr -> Type () == VAR) ||
488 (vr -> Type () == AUX)) &&
489 (vr -> Original () != p -> Var (indVar = vr -> Index ()))) {
490
491 expression *trash = vr;
492 row -> first = p -> Var (indVar);
493 delete trash;
494 }
495
496 // substitute each variable of this row with its newest version
497
498 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
499
500 exprVar * &vc = col -> first;
501 int indVar;
502
503 // substitute variable representing this row with its newest version
504
505 if (((vc -> Type () == VAR) ||
506 (vc -> Type () == AUX)) &&
507 (vc -> Original () != p -> Var (indVar = vc -> Index ()))) {
508
509 expression *trash = vc;
510 col -> first = p -> Var (indVar);
511 delete trash;
512 }
513 }
514 }
515 }
516
517
518 /// compute $y^{lv}$ and $y^{uv}$ for Violation Transfer algorithm
closestFeasible(expression * varind,expression * vardep,CouNumber & left,CouNumber & right) const519 void exprQuad::closestFeasible (expression *varind,
520 expression *vardep,
521 CouNumber &left,
522 CouNumber &right) const {
523
524 fprintf (stderr, "exprQuad::closestFeasible() not available for VT\n");
525 exit (-1);
526 }
527
528
529 /// return l-2 norm of gradient at given point
gradientNorm(const double * x)530 CouNumber exprQuad::gradientNorm (const double *x) {
531
532 CouNumber grad = 0.;
533
534 for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
535
536 CouNumber gradEl = 0.;
537 for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
538 gradEl += col -> second * x [col -> first -> Index ()];
539
540 grad += gradEl * gradEl;
541 }
542
543 return sqrt (grad);
544 }
545
546 /// Simplify expression
simplify()547 expression *exprQuad::simplify () {
548 exprOp::simplify ();
549 return NULL;
550 }
551
552