source: stable/0.2/Couenne/src/expression/operators/exprQuad.cpp @ 159

Last change on this file since 159 was 159, checked in by pbelotti, 11 years ago

created new stable branch 0.2 from trunk (rev. 157)

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