source: trunk/Couenne/src/expression/operators/exprQuad.cpp @ 75

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

fixed segfault with constant objective. Other minor fixes

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