/home/coin/SVN-release/OS-2.2.0/Couenne/src/standardize/decomposeTerm.cpp

Go to the documentation of this file.
00001 /* $Id: decomposeTerm.cpp 230 2009-07-18 11:42:59Z pbelotti $
00002  *
00003  * Name:    decomposeTerm.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: decompose sums and products
00006  *
00007  * (C) Carnegie-Mellon University, 2007. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CouenneProblem.hpp"
00012 #include "exprMul.hpp"
00013 #include "exprPow.hpp"
00014 #include "exprQuad.hpp"
00015 #include "lqelems.hpp"
00016 
00017 //#define DEBUG
00018 
00029 
00030 void CouenneProblem::decomposeTerm (expression *term,
00031                                     CouNumber initCoe,
00032                                     CouNumber &c0,
00033                                     LinMap  &lmap,
00034                                     QuadMap &qmap) {
00035 
00036   switch (term -> code ()) {
00037 
00038     // easy cases ////////////////////////////////////////////////////////////////////////
00039 
00040   case COU_EXPRCONST: 
00041     c0 += initCoe * term -> Value ();
00042     break;
00043 
00044   case COU_EXPRVAR:   
00045     lmap.insert (term -> Index (), initCoe);
00046     break;
00047 
00048   case COU_EXPROPP:   
00049     decomposeTerm (term -> Argument (), -initCoe, c0, lmap, qmap);
00050     break;
00051 
00052   case COU_EXPRSUB:   
00053     decomposeTerm (term -> ArgList () [0],  initCoe, c0, lmap, qmap);
00054     decomposeTerm (term -> ArgList () [1], -initCoe, c0, lmap, qmap);
00055     break;
00056 
00057   case COU_EXPRQUAD: { 
00058 
00059     exprQuad *t = dynamic_cast <exprQuad *> (term -> isaCopy () ? 
00060                                              term -> Copy () : 
00061                                              term);
00062     exprQuad::sparseQ &M = t -> getQ ();
00063 
00064     for (exprQuad::sparseQ::iterator row = M.begin (); 
00065          row != M.end (); ++row) {
00066 
00067       int xind = row -> first -> Index ();
00068 
00069       for (exprQuad::sparseQcol::iterator col = row -> second.begin (); 
00070            col != row -> second.end (); ++col) {
00071         qmap.insert (xind, col -> first -> Index (), initCoe * col -> second);
00072       }
00073     }
00074   } // NO break here, exprQuad generalizes exprGroup
00075 
00076   case COU_EXPRGROUP: { 
00077 
00078     exprGroup *t = dynamic_cast <exprGroup *> (term -> isaCopy () ? 
00079                                                term -> Copy () : 
00080                                                term);
00081     exprGroup::lincoeff &lcoe = t -> lcoeff ();
00082 
00083     //  for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00084     for (int n = lcoe.size (), i=0; n--; i++)
00085       lmap.insert (lcoe [i].first -> Index (), initCoe * lcoe [i].second);
00086 
00087     c0 += initCoe * t -> getc0 ();
00088   } // NO break here, exprGroup generalizes exprSum
00089 
00090   case COU_EXPRSUM: { 
00091 
00092     expression **al = term -> ArgList ();
00093     for (int i = term -> nArgs (); i--;)
00094       decomposeTerm (*al++, initCoe, c0, lmap, qmap);
00095 
00096   } break;
00097 
00098   // not-so-easy cases /////////////////////////////////////////////////////////////////
00099   //
00100   // cannot add terms as it may fill up the triplet
00101 
00102   case COU_EXPRMUL: { 
00103 
00104     std::map <int, CouNumber> indices;
00105     CouNumber coe = initCoe;
00106 
00107     // return list of variables (some of which auxiliary)
00108     flattenMul (term, coe, indices);
00109 
00110 #ifdef DEBUG
00111     printf ("from flattenmul: [%g] ", coe);
00112     for (std::map <int, CouNumber>::iterator itt = indices.begin ();
00113          itt != indices.end(); ++itt)
00114       printf (" %d,%g",
00115               itt -> first,
00116               itt -> second);
00117     printf ("\n");
00118 #endif
00119 
00120     // based on number of factors, decide what to return
00121     switch (indices.size ()) {
00122 
00123     case 0: // no variables in multiplication (hmmm...)
00124       c0 += coe;
00125       break;
00126 
00127     case 1: { // only one term (may be with !=1 exponent)
00128 
00129       std::map <int, CouNumber>::iterator one = indices.begin ();
00130       int       index = one -> first;
00131       CouNumber expon = one -> second;
00132 
00133       if      (fabs (expon - 1) < COUENNE_EPS) lmap.insert (index, coe);
00134       else if (fabs (expon - 2) < COUENNE_EPS) qmap.insert (index, index, coe);
00135       else {
00136         exprAux *aux = addAuxiliary 
00137           (new exprPow (new exprClone (Var (index)),
00138                         new exprConst (expon)));
00139 
00140         //linsert (lmap, aux -> Index (), initCoe); // which of these three is correct?
00141         //linsert (lmap, aux -> Index (), initCoe * coe);
00142         lmap.insert (aux -> Index (), coe);
00143       }
00144     } break;
00145 
00146     case 2: { // two terms
00147 
00148       int ind0, ind1;
00149 
00150       std::map <int, CouNumber>::iterator one = indices.begin (), 
00151         two = one;
00152       ++two; // now "two" points to the other variable
00153 
00154       // first variable
00155       if (fabs (one -> second - 1) > COUENNE_EPS) {
00156         exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (one -> first)),
00157                                                   new exprConst (one -> second)));
00158         ind0 = aux -> Index ();
00159       } else ind0 = one -> first;
00160 
00161       // second variable
00162       if (fabs (two -> second - 1) > COUENNE_EPS) {
00163         exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (two -> first)),
00164                                                   new exprConst (two -> second)));
00165         ind1 = aux -> Index ();
00166       } else ind1 = two -> first;
00167 
00168       qmap.insert (ind0, ind1, coe);
00169     } break;
00170 
00171     default: { 
00172 
00173       // create new auxiliary variable containing product of 3+ factors
00174 
00175       expression **al = new expression * [indices.size ()];
00176       std::map <int, CouNumber>::iterator one = indices.begin ();
00177 
00178       for (int i=0; one != indices.end (); ++one, i++) 
00179         if (fabs (one -> second - 1) > COUENNE_EPS) {
00180           exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (one -> first)),
00181                                                     new exprConst (one -> second)));
00182           al [i] = new exprClone (aux);
00183         } else al [i] = new exprClone (Var (one -> first));
00184 
00185       // TODO: when we have a convexification for \prod_{i \in I}...
00186       //      exprAux *aux = addAuxiliary (new exprMul (al, indices.size ()));
00187 
00188       exprMul *mul = new exprMul (al, indices.size ());
00189       exprAux *aux = mul -> standardize (this);
00190       lmap.insert (aux -> Index (), coe);
00191 
00192     } break;
00193     }
00194 
00195   } break; // end of case COU_EXPRMUL
00196 
00197   case COU_EXPRPOW: { // expression = f(x)^g(x) ////////////////////////////////////////////////
00198 
00199     expression **al = term -> ArgList (); 
00200 
00201     if (al [1] -> Type () != CONST) { 
00202 
00203       // non-constant exponent, standardize the whole term and add
00204       // linear component (single aux)
00205 
00206       expression *aux = term -> standardize (this);
00207       if (!aux) aux = term;
00208       lmap.insert (aux -> Index (), initCoe);
00209 
00210     } else { // this is of the form f(x)^k.  If k=2, return square. If
00211              // k=1, return var. Otherwise, generate new auxiliary.
00212 
00213       expression *aux = (*al) -> standardize (this);
00214 
00215       if (!aux)
00216         aux = *al; // it was a simple variable, and was not standardized.
00217 
00218       CouNumber expon = al [1] -> Value ();
00219       int ind = aux -> Index ();
00220 
00221       if      (fabs (expon - 1) < COUENNE_EPS) lmap.insert (ind, initCoe);
00222       else if (fabs (expon - 2) < COUENNE_EPS) qmap.insert (ind, ind, initCoe);
00223       else {
00224         exprAux *aux2 = addAuxiliary 
00225           (new exprPow (new exprClone (aux), new exprConst (expon))); // TODO: FIX!
00226         lmap.insert (aux2 -> Index (), initCoe);
00227       }
00228     }
00229   } break;
00230 
00231   default: { 
00232 
00233     expression *aux = term -> standardize (this);
00234     if (!aux) 
00235       aux = term;
00236     lmap.insert (aux -> Index (), initCoe);
00237   } break;
00238   }
00239 }

Generated on Thu Aug 5 03:02:57 2010 by  doxygen 1.4.7