/home/coin/SVN-release/OS-2.4.1/Couenne/src/expression/operators/exprQuad.cpp

Go to the documentation of this file.
00001 /* $Id: exprQuad.cpp 597 2011-06-02 10:09:33Z pbelotti $
00002  *
00003  * Name:    exprQuad.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: implementation of some methods for exprQuad
00006  *
00007  * (C) Carnegie-Mellon University, 2006-08.
00008  * This file is licensed under the Eclipse Public License (EPL)
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 //#define DEBUG
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       //coe /= 2.;
00062 
00063       // pick smaller index as row reference
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 ()) { // add new row
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 { // insert element into row
00084 
00085       rowMap::iterator cell = rowp -> second.find (varJ);
00086 
00087       if (cell == rowp -> second.end ()) { // normal case, add entry
00088 
00089         std::pair <exprVar *, CouNumber> newcell (varJ, coe);
00090         rowp -> second.insert (newcell);
00091 
00092       } else { // strange, but add coefficient
00093 
00094         if (fabs (cell -> second += coe) < COUENNE_EPS)
00095           // eliminate element of map if null coefficient
00096           rowp -> second.erase (cell); 
00097       }
00098     }
00099   }
00100 
00101   // transform maps into vectors
00102 
00103   for (matrixMap::iterator row = qMap.begin (); row != qMap.end (); ++row) {
00104 
00105     sparseQcol line;
00106 
00107     // insert first element in bound map
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       // insert second element in bound map
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                         //(dynamic_cast <exprVar *> (i -> first -> clone (d)), i -> second));
00147                         (new exprVar (i -> first -> Index (), d), i -> second));
00148 
00149     matrix_.push_back (std::pair <exprVar *, sparseQcol> 
00150                        //dynamic_cast <exprVar *> (row -> first -> clone (d)), column));
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   //if (code () == COU_EXPRQUAD)
00181   if (matrix_.size () > 0)
00182     out << '(';
00183 
00184   // print linear and nonquadratic part
00185   exprGroup::print (out, descend);
00186 
00187   int noperands = 0;
00188 
00189   for (size_t n = matrix_.size (), i=0; n--; i++) {
00190     //sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
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       //sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
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   //if (code () == COU_EXPRGROUP)
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   // derive linear part (obtain constant)
00234   for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00235     c0 += el -> second;
00236 
00237   // derive quadratic part (obtain linear part)
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   // derive nonlinear sum
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   // special cases
00276 
00277   // 1) no linear part
00278   if (lmap.empty ()) {
00279 
00280     // and no nonlinear part either
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     //    if (matrix_.size() > e.matrix_.size()) return  1;
00326     //    int xind = row -> first -> Index ();
00327     //    CouNumber x = (*(row -> first)) ();
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  // coefficient fractional, check all is fixed and product is integer
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       //if (!intI && !intJ &&  intCoe) ; // check x y fixed int
00435       //if (!intI &&  intJ &&  intCoe) ; // check x   fixed int
00436       //if ( intI && !intJ &&  intCoe) ; // check   y fixed int
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       //fprintf (stderr, "Didn't fix exprQuad::replace() yet");
00459       //exit (-1);
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         //fprintf (stderr, "Didn't fix exprQuad::replace() yet");
00469         //exit (-1);
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     // substitute variable representing this row with its newest version
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     // substitute each variable of this row with its newest version
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       // substitute variable representing this row with its newest version
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 

Generated on Thu Nov 10 03:05:44 2011 by  doxygen 1.4.7