/home/coin/SVN-release/OS-2.0.1/Couenne/src/standardize/splitAux.cpp

Go to the documentation of this file.
00001 /* $Id: splitAux.cpp 154 2009-06-16 18:52:53Z pbelotti $ */
00002 /*
00003  * Name:    splitAux.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: extract auxiliary variable from implicit equality constraint
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 
00013 #include "CoinHelperFunctions.hpp"
00014 
00015 #include "exprClone.hpp"
00016 #include "exprSum.hpp"
00017 #include "exprMul.hpp"
00018 #include "exprGroup.hpp"
00019 #include "exprQuad.hpp"
00020 #include "lqelems.hpp"
00021 
00022 //#define DEBUG
00023 
00027 void elementBreak (expression *, int &, CouNumber &);
00028 
00029 
00032 
00033 int CouenneProblem::splitAux (CouNumber rhs, expression *body, expression *&rest, bool *wentAux) {
00034 
00035   int auxInd = -1,          // index of the auxiliary to be extracted
00036     code = body -> code (); // type of expression
00037 
00038   expression **alist = body -> ArgList ();
00039 
00040 #ifdef DEBUG
00041   printf ("|||||||||| Splitting "); body -> print (); printf ("\n");
00042 #endif
00043 
00044   switch (code) { // constraint h(x) = 0 may be in different forms:
00045                   // subtraction, sum, linear group
00046 
00048 
00049   case COU_EXPRSUB: {
00050 
00051     // simplest case, we have w-f(x)=rhs or f(x)-w=rhs or f(x)-g(x)=rhs
00052 
00053     int pos = 0;
00054     CouNumber coeff = 1;
00055 
00056     auxInd = (*alist) -> Index ();
00057 
00058     if (auxInd < 0)
00059       elementBreak (*alist, auxInd, coeff); // check first element 
00060 
00061     if ((auxInd < 0) || 
00062         wentAux [auxInd] || 
00063         (alist [1] -> dependsOn (auxInd, TAG_AND_RECURSIVE) >= 1)) {
00064 
00065       auxInd = (alist [1]) -> Index ();
00066 
00067       if (auxInd < 0)
00068         elementBreak (alist [1], auxInd, coeff); // check second element
00069       else coeff = 1;
00070 
00071       if ((auxInd < 0) ||       // no index found
00072           wentAux [auxInd] ||   // or, found but invalid
00073           (alist [0] -> dependsOn (auxInd, TAG_AND_RECURSIVE) >= 1)) 
00074                                 // or, variable depends upon itself
00075         return -1;
00076 
00077       pos = 1;
00078     }
00079 
00081 
00082     // what remains is the "independent" expression
00083     expression *clone = alist [1 - pos] -> clone (&domain_); 
00084 
00085     expression *auxdef = 
00086       (fabs (coeff - 1) < COUENNE_EPS) ?          // if coefficient is 1
00087       (clone) :                                   //     do nothing
00088       new exprMul (new exprConst (coeff), clone); //     else multiply it by coefficient
00089 
00090     rest = (fabs (rhs) < COUENNE_EPS) ?                              // no extra constant?
00091       (auxdef) :                                                     // just put other argument
00092       (new exprSum (auxdef, new exprConst ((pos==1) ? -rhs : rhs))); // otherwise sum it with \pm rhs
00093 
00094   } break;
00095 
00097 
00098   case COU_EXPRQUAD:
00099   case COU_EXPRGROUP:
00100   case COU_EXPRSUM: {
00101 
00102     // an expr{Sum,Group,Quad}. Decompose the expression and
00103     // re-assemble (to look for right variable)
00104 
00105     // data structure to be used below if there is a linear term.
00106     // which specifies position within arrays (negative for linear
00107     // part of exprGroup, positive for all elements of exprSum)
00108 
00109     int maxindex = -1, nlin = 0, which = 1;
00110     CouNumber c0 = 0., auxcoe = 1;
00111     bool which_was_set = false;
00112 
00113     // check indices of linear part /////////////////////////////////
00114 
00115     if (code != COU_EXPRSUM) {
00116 
00117       exprGroup *egBody = dynamic_cast <exprGroup *> (body -> isaCopy () ? 
00118                                                       body -> Copy () : 
00119                                                       body);
00120       exprGroup::lincoeff &lcoe = egBody -> lcoeff ();
00121 
00122       // import exprGroup linear structure
00123 
00124       c0 = egBody -> getc0 ();
00125 
00126       for (int i=0, n = lcoe.size (); n--; i++, nlin++) {
00127 
00128         int j = lcoe [i]. first -> Index ();
00129 
00130         //lincoe [i] = lcoe [i]. second;
00131 
00132         // prefer non-integer. If integer, only take it if none chosen yet
00133         if ((!(lcoe [i].first -> isInteger ()) || (which==1)) &&  
00134             (j > maxindex) &&
00135             !(wentAux [j]) && 
00136             (fabs (lcoe [i]. second) > COUENNE_EPS)) {
00137 
00138           // fake cut in linind and check dependence. Only mark if
00139           // dependsOn() gives 0
00140 
00141           exprVar *saveVar = lcoe [i].first;
00142           //lcoe [i].first = new exprVar (nVars ()); // better use index -1
00143           lcoe [i].first = new exprVar (-1);
00144 
00145           if (body -> dependsOn (j, TAG_AND_RECURSIVE) == 0) {
00146 
00147 #ifdef DEBUG
00148             printf ("body does not depend on x%d: ", j);
00149             body -> print ();
00150             printf ("\n");
00151 #endif
00152             // mark which with negative number
00153             which    = - nlin - 1;
00154             auxcoe   = lcoe [i]. second;
00155             maxindex = j;
00156           }
00157 
00158           delete lcoe [i].first;
00159           lcoe [i].first = saveVar;
00160         }
00161       }
00162     }
00163 
00164     if (which != 1) 
00165       which_was_set = true;
00166     else which = -1;
00167 
00168     // check indices of elements of (nonlinear) sum /////////////////////////////////
00169 
00170     for (int i = body -> nArgs (); i--;) {
00171 
00172       CouNumber coeff = 1;
00173       int index = alist [i] -> Index ();
00174 
00175       if (index < 0)
00176         elementBreak (alist [i], index, coeff);
00177 
00178       if ((index > maxindex) &&
00179           !(wentAux [index]) &&
00180           (fabs (coeff) > COUENNE_EPS)) {
00181 
00182         // fake a cut in the arglist and check
00183 
00184         expression *cut = alist [i];
00185         alist [i] = new exprConst (0.);
00186 
00187         // not enough... check now linear (and quadratic!) terms 
00188 
00189         if (body -> dependsOn (index, TAG_AND_RECURSIVE) == 0) {
00190 
00191           maxindex = index;
00192           which    = i;
00193           auxcoe   = coeff;
00194         }
00195 
00196         delete alist [i];
00197         alist [i] = cut;
00198       }
00199     }
00200 
00201     if (!which_was_set && (which == -1)) // which has not been set
00202       return -1;
00203 
00205 
00206     if (maxindex < 0) break; // no substitute found ==> no hidden auxiliary
00207 
00208     // create a new exprGroup, exprQuad, or exprSum with all elements but the
00209     // extracted auxiliary
00210 
00211     // start with exprSum
00212     int nargs = body -> nArgs ();
00213     expression **newarglist;
00214 
00215 #ifdef DEBUG
00216     printf (" [[ind %d, coe %.1g, wh %d, nargs %d, nlin %d]] ", 
00217             maxindex, auxcoe, which, nargs, nlin);
00218 #endif
00219 
00220     if (nargs > 0) { // there is an element in the nonlinear sum to be drawn
00221 
00222       int j, mid = (which < 0) ? nargs : which;
00223 
00224       newarglist = new expression * [nargs + 1];
00225 
00226       for (j=0; j<mid;   j++) newarglist [j]   = alist [j] -> clone (&domain_);
00227       for (j++; j<nargs; j++) newarglist [j-1] = alist [j] -> clone (&domain_);
00228 
00229       // nl arglist is done, later decide whether to incorporate it as
00230       // it is or with a coefficient
00231 
00232     } else { // no nonlinear arguments, or the only one is the new aux
00233 
00234       nargs++; // !!!!!!!!!!!!!!!!!!!!!!!!!
00235       newarglist  = new expression *;
00236       *newarglist = new exprConst (0.);
00237     }
00238 
00239     // form rhs linear part ////////////////////////////////////////////////////
00240 
00241     int       *linind2 = NULL;
00242     CouNumber *lincoe2 = NULL;
00243 
00244     // in case this was (and will be) an exprQuad
00245     int *qindI = NULL, 
00246         *qindJ = NULL;
00247     CouNumber *qcoe = NULL;
00248 
00249     if (nlin > 0) { // there is an element in the linear sum to be drawn
00250 
00251       exprGroup *egBody = dynamic_cast <exprGroup *> (body -> isaCopy () ? 
00252                                                       body -> Copy () : 
00253                                                       body);
00254       exprGroup::lincoeff &lcoe = egBody -> lcoeff ();
00255 
00256       int mid = (which >= 0) ? nlin : - which - 1;
00257 
00258       linind2 = new int       [nlin + 1];
00259       lincoe2 = new CouNumber [nlin + 1];
00260 
00261       register int j;
00262 
00263 #ifdef DEBUG
00264       //for (j=0; j<mid;  j++) printf ("{%g x%d} ", lincoe [j], linind [j]);
00265       //for (j++; j<nlin; j++) printf ("{%g x%d} ", lincoe [j], linind [j]);
00266 #endif
00267 
00268       CouNumber divider = -1. / auxcoe;
00269 
00270       for (j=0; j<mid;  j++){linind2[j]  =lcoe[j].first->Index();lincoe2[j]  =divider*lcoe[j].second;}
00271       for (j++; j<nlin; j++){linind2[j-1]=lcoe[j].first->Index();lincoe2[j-1]=divider*lcoe[j].second;}
00272 
00273       linind2 [j-1] = -1; // terminate list of indices
00274 
00275 #ifdef DEBUG
00276       for (j=0; j<mid;  j++) printf ("<%g x%d> ", lincoe2 [j],   linind2 [j]);
00277       for (j++; j<nlin; j++) printf ("<%g x%d> ", lincoe2 [j-1], linind2 [j-1]);
00278 #endif
00279 
00280       if (code == COU_EXPRQUAD) { // copy quadratic elements
00281 
00282         exprQuad *eq = dynamic_cast <exprQuad *> (body -> isaCopy () ? 
00283                                                   body -> Copy () : 
00284                                                   body);
00285 
00286         int nqt = eq -> getnQTerms (), j=0;
00287 
00288         qindI = new int [nqt];
00289         qindJ = new int [nqt];
00290         qcoe  = new CouNumber [nqt];
00291 
00292         exprQuad::sparseQ &M = eq -> getQ ();
00293 
00294         for (exprQuad::sparseQ::iterator row = M.begin (); 
00295              row != M.end (); ++row) {
00296 
00297           int xind = row -> first -> Index ();
00298 
00299           for (exprQuad::sparseQcol::iterator col = row -> second.begin (); 
00300                col != row -> second.end (); ++col, ++j) {
00301 
00302             qindI [j] = xind;
00303             qindJ [j] = col -> first -> Index ();
00304             qcoe  [j] = divider * col -> second;
00305           }
00306         }
00307       }
00308 
00309       // nl arglist is done, later decide whether to incorporate it as
00310       // it is or with a coefficient
00311     }
00312 
00313     // the extracted index is one term of...
00314     if (which >= 0) --nargs; // ...the nonlinear sum
00315     else            --nlin;  // ...the linear part
00316 
00317 #ifdef DEBUG
00318     printf ("\n::: auxcoe %g, rhs %g, lin %d, nl %d\n", auxcoe, rhs, nlin, nargs);
00319 #endif
00320 
00321     // all is ready to take the independent stuff to the other side of
00322     // the inequality.
00323 
00324     if ((code == COU_EXPRQUAD) || 
00325         (code == COU_EXPRGROUP) && (nlin > 0)) { 
00326 
00327       // an exprGroup with at least one linear term left
00328       //
00329       // build new vectors for index and coeff. Two cases:
00330       //
00331       // 1)  f(x) + c0 -  w = rhs   =====>   w =       f(x) + c0 - rhs
00332       // 2)  f(x) + c0 + aw = rhs   =====>   w = -1/a (f(x) + c0 - rhs), a != -1
00333 
00334       if (fabs (auxcoe + 1) < COUENNE_EPS) {
00335 
00336         //std::vector <std::pair <exprVar *, CouNumber> >
00337         exprGroup::lincoeff lcoeff;
00338         indcoe2vector (linind2, lincoe2, lcoeff);
00339 
00340         if (code == COU_EXPRGROUP)
00341           rest = new exprGroup (c0-rhs, lcoeff, newarglist, nargs);
00342 
00343         else {
00344           std::vector <quadElem> qcoeff;
00345           indcoe2vector (qindI, qindJ, qcoe, qcoeff);
00346           rest = new exprQuad  (c0-rhs, lcoeff, qcoeff, newarglist, nargs);
00347         }
00348       }
00349       else {
00350 
00351         expression **mullist = new expression * [1];
00352 
00353         // only nl term is constant
00354         if ((nargs <= 1) && ((*newarglist) -> Linearity () <= CONSTANT)) {
00355           *mullist = new exprConst ((*newarglist) -> Value ());
00356           //delete *newarglist;
00357           delete [] newarglist;
00358         }
00359         else if ((nargs <= 1) && 
00360                  ((*newarglist) -> code () == COU_EXPROPP) &&
00361                  (fabs (auxcoe - 1) < COUENNE_EPS))
00362           //*mullist = new exprSum (newarglist, nargs);
00363           *mullist = (*newarglist) -> Argument ();//, nargs);
00364         else { // multiple nonlinear terms, multiply them by -1/auxcoe
00365 
00366           if (nargs <= 1)
00367             *mullist = new exprMul (new exprConst (-1. / auxcoe), 
00368                                     *newarglist); // !!! one more leak?
00369           else 
00370             *mullist = new exprMul (new exprConst (-1. / auxcoe), 
00371                                     new exprSum (newarglist, nargs));
00372         }
00373 
00374         std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
00375         indcoe2vector (linind2, lincoe2, lcoeff);
00376 
00377         // final outcome: -1/a (f(x) + c0 - rhs)
00378         if (code == COU_EXPRGROUP)
00379           rest = new exprGroup ((rhs - c0) / auxcoe, lcoeff,        mullist,1);
00380 
00381         else {
00382           std::vector <quadElem> qcoeff;
00383           indcoe2vector (qindI, qindJ, qcoe, qcoeff);
00384           rest = new exprQuad  ((rhs - c0) / auxcoe, lcoeff, qcoeff, mullist,1);
00385         }
00386       }
00387     }
00388     else { // simple exprSum
00389 
00390       if (fabs (rhs) > COUENNE_EPS) { // have to add constant to exprSum
00391 
00392         if ((nargs == 1) && ((*newarglist) -> Type () == CONST)) {
00393 
00394           CouNumber val = (*newarglist) -> Value () - rhs;
00395           delete *newarglist;
00396           *newarglist = new exprConst (val);
00397         } 
00398         else newarglist [nargs++] = new exprConst (-rhs);
00399       }
00400 
00401       // now exprSum is complete with -rhs. Send it to right hand side
00402 
00403       expression *auxDef;
00404 
00405       if (nargs==1) {
00406         auxDef = *newarglist;
00407         delete [] newarglist;
00408       } else auxDef = new exprSum (newarglist, nargs);
00409 
00410       if (fabs (auxcoe + 1) < COUENNE_EPS)
00411         rest = auxDef;
00412       else if ((fabs (auxcoe - 1) < COUENNE_EPS) && 
00413                (auxDef -> code () == COU_EXPROPP))
00414         rest = auxDef -> Argument ();
00415       else // TODO: check if auxdef is an exprOpp or an exprMul
00416            // (k*f(x)) and -1/auxcoe simplifies
00417         rest = new exprMul (new exprConst (-1./auxcoe), 
00418                             ((auxDef -> Type () == AUX) ||
00419                              (auxDef -> Type () == VAR)) ? new exprClone (auxDef) : auxDef);
00420     }
00421 
00422     if (linind2) {
00423       delete [] linind2;
00424       delete [] lincoe2;
00425     }
00426 
00427 #ifdef DEBUG
00428     printf ("gotten "); rest -> print (); printf ("\n");
00429 #endif
00430 
00431     auxInd = maxindex;
00432 
00433   } break;
00434 
00435   default: break;
00436   } // end switch () ////////////////////////////////////////////////////////
00437 
00438   if ((auxInd < 0) || (wentAux [auxInd]))
00439     return -1;
00440 
00441   // we have a variable, meaning the constraint body is of the form
00442   // w +/- f(x) or f(x) +/- w
00443 
00444   // standardize remaining of the expression
00445 
00446 #ifdef DEBUG
00447   printf ("standardize rest (2nd level) "); fflush (stdout);
00448   rest -> print (); printf (" [body = ");
00449   body -> print (); printf ("]");
00450 #endif
00451 
00452   // second argument is false to tell standardize not to create a new
00453   // auxiliary variable (we are creating it ourselves, it's aux)
00454   exprAux *aux = rest -> standardize (this, false);
00455 
00456 #ifdef DEBUG
00457   printf (" {body = "); body -> print (); printf ("} ");
00458   if (aux) {
00459     printf (" --- aux = "); 
00460     aux -> print ();
00461     printf (" := ");
00462     aux -> Image () -> print ();
00463   }
00464 #endif
00465 
00466   if (aux) {
00467     //delete rest;
00468     //rest = aux -> Image () -> clone (&domain_);
00469     rest = aux -> Image ();// -> clone (&domain_);
00470     aux -> Image (NULL);
00471     delete aux;
00472   }
00473 
00474 #ifdef DEBUG
00475   printf (" ==> "); rest -> print ();
00476   printf (" and "); fflush (stdout);
00477   rest -> print (); printf ("\n");
00478 #endif
00479 
00480   return auxInd;
00481 }

Generated on Thu Oct 8 03:02:57 2009 by  doxygen 1.4.7