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

Go to the documentation of this file.
00001 /* $Id: standardize.cpp 752 2011-08-08 03:45:07Z pbelotti $
00002  *
00003  * Name:    standardize.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: standardize all expressions in a problem
00006  *
00007  * (C) Carnegie-Mellon University, 2006-11.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include <vector>
00012 
00013 #include "CoinHelperFunctions.hpp"
00014 
00015 #include "BonBabSetupBase.hpp"
00016 
00017 #include "CouenneTypes.hpp"
00018 #include "CouenneExpression.hpp"
00019 #include "CouenneExprIVar.hpp"
00020 #include "CouenneExprSub.hpp"
00021 #include "CouenneExprClone.hpp"
00022 #include "CouenneProblem.hpp"
00023 #include "CouenneProblemElem.hpp"
00024 #include "CouenneDepGraph.hpp"
00025 
00026 using namespace Ipopt;
00027 using namespace Couenne;
00028 
00029 void replace (CouenneProblem *p, int wind, int xind);
00030 
00032 
00033 bool CouenneProblem::standardize () {
00034 
00035   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00036     printf ("Reformulation. current point: %d vars -------------------\n", 
00037             domain_.current () -> Dimension ());
00038     for (int i=0; i<domain_.current () -> Dimension (); i++)
00039       printf ("%3d %20g [%20g %20g]\n", i, domain_.x (i), domain_.lb (i), domain_.ub (i));
00040   }
00041 
00042   bool retval = true;
00043 
00044   // create dependence graph to assign an order to the evaluation (and
00045   // bound propagation, and -- in reverse direction -- implied bounds)
00046   graph_ = new DepGraph;
00047 
00048   for (std::vector <exprVar *>::iterator i = variables_ . begin ();
00049        i != variables_ . end (); ++i)
00050     graph_ -> insert (*i);
00051 
00052   // allocate space in auxiliaries_ from commonexprs_
00053 
00054   //int initVar = variables_ . size () - commonexprs_ . size ();
00055 
00056   //std::set <int> DefVarDiffSet;
00057 
00058   // // DEFINED VARIABLES -----------------------------------------------------------------------
00059 
00060   // // standardize initial aux variables (aka defined variables, aka
00061   // // common expression)
00062 
00063   // if (commonexprs_.size ()) jnlst_ -> Printf (J_ALL, J_REFORMULATE,
00064   //                                          "%d common exprs, initVar = %d = %d - %d\n", 
00065   //                                          commonexprs_ . size (), 
00066   //                                          initVar, 
00067   //                                          variables_   . size (), 
00068   //                                          commonexprs_ . size ());
00069 
00070   // for (std::vector <expression *>::iterator i = commonexprs_ . begin ();
00071   //      i != commonexprs_ . end (); ++i) {
00072 
00073   //   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00074   //     printf ("\n=====> [1] stdz common expr [%d] := ", initVar); fflush (stdout);
00075   //     (*i) -> print (); printf ("\n");
00076   //   }
00077 
00078   //   exprAux    *naux = (*i) -> standardize (this, false);
00079   //   expression *img  = (*i);
00080 
00081   //   if (naux)
00082   //     img = naux -> Image ();
00083 
00084   //   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00085   //     printf ("\n=====> [2] "); fflush (stdout);
00086   //     if (*i)   (*i) -> print (); else printf ("(null)"); printf (" [*i] ");   fflush (stdout);
00087   //     if (naux) naux -> print (); else printf ("(null)"); printf (" [naux] "); fflush (stdout);
00088   //     if (img)  img  -> print (); else printf ("(null)"); printf (" [img]\n");
00089   //   }
00090 
00091   //   // trick to obtain same index as common expression: create exprAux
00092   //   // with index initVar and replace occurrences with address of
00093   //   // newly created exprAux through auxiliarize()
00094 
00095   //   exprVar *newvar;
00096   //   //img -> isInteger () ? exprAux::Integer : exprAux::Continuous);
00097 
00098   //   //auxiliarize (newvar); // puts newvar at right position in variables_
00099 
00100   //   if (((*i) -> Type () == VAR) ||
00101   //    ((*i) -> Type () == AUX)) {
00102 
00103   //     newvar = new exprAux (img, initVar, 1 + img -> rank (), exprAux::Unset, &domain_);
00104   //     replace (this, initVar, img -> Index ());
00105   //     //auxiliarize (variables_ [initVar], 
00106   //     //variables_ [img -> Index ()]);
00107 
00108   //     //delete variables_ [initVar];
00109 
00110   //     variables_ [initVar] = newvar;
00111   //     variables_ [initVar] -> zeroMult ();
00112 
00113   //   } else {
00114 
00115   //     if (img -> dependsOn (&initVar, 1, TAG_AND_RECURSIVE)) {
00116 
00117   //    //printf ("depends! "); img -> print (); 
00118 
00119   //    expression *diff = new exprSub (new exprClone (variables_ [initVar]), img);
00120 
00121   //    //printf ("; diff = "); diff -> print ();
00122 
00123   //    exprAux *diffAux = diff -> standardize (this, false);
00124 
00125   //    //printf ("; aux: "); if (diffAux) diffAux -> print (); 
00126 
00127   //    //if (diffAux)
00128   //    exprAux *newAux = addAuxiliary (diff);
00129 
00130   //    //printf ("; real aux: "); if (newAux) newAux -> print (); putchar ('\n');
00131         
00132   //    //Lb (newAux -> Index ()) = 
00133   //    //Ub (newAux -> Index ()) = 0.;
00134 
00135   //    DefVarDiffSet. insert (newAux -> Index ());
00136 
00137   //     } else {
00138 
00139   //    newvar = new exprAux (img, initVar, 1 + img -> rank (), exprAux::Unset, &domain_);
00140   //    //replace (this, initVar, newvar -> Index ());
00141 
00142   //    auxiliarize (newvar);
00143 
00144   //    //delete variables_ [initVar];
00145   //    variables_ [initVar] = newvar;
00146 
00147   //    graph_ -> insert (newvar);
00148   //    //if (naux) 
00149   //    graph_ -> erase (naux);
00150   //     }
00151   //   }
00152 
00153   //   //variables_ . erase (variables_ . end () - 1);
00154 
00155   //   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00156   //     if (naux) {
00157   //    printf ("done: "); fflush (stdout);
00158   //    naux -> print ();
00159   //    printf (" := "); fflush (stdout);
00160   //    naux -> Image () -> print (); printf ("\n..."); fflush (stdout);
00161   //     } else if (*i) {
00162   //    (*i) -> print ();
00163   //    //printf (" := "); fflush (stdout);
00164   //    //aux -> Image () -> print (); 
00165   //    printf ("\n");
00166   //     } else printf ("[n]aux NULL!\n");
00167   //   }
00168 
00169   //   initVar++;
00170   // }
00171 
00172   // OBJECTIVES ------------------------------------------------------------------------------
00173 
00174   for (std::vector <CouenneObjective *>::iterator i = objectives_.begin ();
00175        i != objectives_.end (); ++i) {
00176 
00177     if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00178       printf ("Objective [code: %d]", (*i) -> Body () -> code ()); 
00179       (*i) -> print ();
00180     }
00181 
00182     exprAux *aux = (*i) -> standardize (this);
00183 
00184     if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00185       printf (" objective "); (*i) -> print (); 
00186       if (aux) {printf (" admits aux "); aux -> print ();}
00187     }
00188 
00189     if (aux) {
00190       //delete ((*i) -> Body ());
00191       (*i) -> Body (new exprClone (aux));
00192     }
00193 
00194     if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00195       printf (". New obj: "); (*i) -> print (); printf ("\n");
00196     }
00197   }
00198 
00199   // Constraints ----------------------------------------------------------------------------
00200 
00201   // commuted_ is an array with a flag for each original variable,
00202   // which is true at position i if initially original variable x_i
00203   // became auxiliary (because of constraint 
00204   // 
00205   // x_{k+1} + f(x_1,x_2...,x_k} <=/=/>= 0
00206   //
00207   // becoming
00208   //
00209   // x_{k+1} <=/=/>= - f(x_1,x_2...,x_k}
00210 
00211   commuted_ = new bool [nVars ()];
00212   CoinFillN (commuted_, nVars (), false);
00213 
00214   std::vector <std::vector <CouenneConstraint *>::iterator> iters2erase;
00215 
00216   for (std::vector <CouenneConstraint *>::iterator i = constraints_.begin (); 
00217        i != constraints_.end (); ++i) {
00218 
00219     if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00220       printf ("\nReformulating constraint: "); 
00221       (*i) -> print ();
00222     }
00223 
00224     CouNumber 
00225       conLb = (*((*i) -> Lb ())) (),
00226       conUb = (*((*i) -> Ub ())) ();
00227 
00228     // sanity check: if (due to bad model or redundancies) the
00229     // constraint's body is constant, delete it -- check if it's
00230     // within bounds.
00231 
00232     expression *eBody = (*i) -> Body ();
00233 
00234     if (eBody -> Linearity () <= CONSTANT) {
00235 
00236       CouNumber bodyVal = (*eBody)();
00237 
00238       if ((bodyVal < conLb - COUENNE_BOUND_PREC) ||
00239           (bodyVal > conUb + COUENNE_BOUND_PREC)) { // all variables eliminated, but out of bounds
00240         
00241         jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, 
00242                           "Constraint: all variables eliminated, but value %g out of bounds [%g,%g]: ",                           
00243                           bodyVal, conLb, conUb);
00244 
00245         if (jnlst_ -> ProduceOutput (J_SUMMARY, J_PROBLEM))
00246           (*i) -> print ();
00247 
00248         retval = false;
00249         break;
00250 
00251       } else {
00252 
00253         iters2erase.push_back (i);
00254         continue; // all variables eliminated and constraint is redundant
00255       }
00256     }
00257 
00258     exprAux *aux = (*i) -> standardize (this);
00259 
00260     if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00261       printf (" reformulated: aux w[%d] ", aux ? (aux -> Index ()) : -2); 
00262       (*i) -> print ();
00263     }
00264 
00265     if (aux) { // save if standardized
00266 
00267       // this is a top level auxiliary, i.e., an auxiliary whose node
00268       // in the DAG stands at the maximum level -- no other variable
00269       // depends on it as it is the lhs of a constraint.
00270 
00271       aux -> top_level () = true;
00272 
00273       // this constraint f(x)<=b is therefore replaced by w<=b, where
00274       // w is the top-level auxiliary associated with f(x)
00275 
00276       //printf ("delete %x: ", ((*i) -> Body ())); ((*i) -> Body ()) -> print ();
00277       //printf ("\n"); 
00278       //delete ((*i) -> Body ());
00279 
00280       (*i) -> Body (new exprClone (aux));
00281 
00282       //      con2.push_back (*i);
00283     }
00284     else {
00285 
00286       // left-hand side not reformulated, therefore this is probably
00287       // an affine expression
00288 
00289       //CouNumber lb, ub;
00290 
00291       //printf ("let's see what's left of the body: "); fflush (stdout);
00292       //(*i) -> Body () -> print (); printf ("\n");
00293 
00294       // (*i) -> Body () -> getBounds (lb, ub);
00295 
00296       // if ((((*((*i) -> Lb ())) ()) > ub + COUENNE_EPS) ||
00297       //          (((*((*i) -> Ub ())) ()) < lb - COUENNE_EPS)) {
00298 
00299       //        jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, "found infeasible constraint [%g,%g] vs [%g,%g]\n", 
00300       //                          lb, ub,
00301       //                          ((*((*i) -> Lb ())) ()),
00302       //                          ((*((*i) -> Ub ())) ()));
00303 
00304       //        if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE))
00305       //          (*i) -> print ();
00306 
00307       //        retval = false;
00308       // }
00309 
00310       // delete constraint as no aux associated with it. It is a
00311       // linear constraint now (after reformulation), so it will be
00312       // inserted as a linearization of the aux's definition (an
00313       // affine function).
00314 
00315       iters2erase.push_back (i);
00316     }
00317 
00318     //(*i) -> Body () -> realign (this);
00319 
00320     if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00321       printf (" --> "); (*i) -> print (); printf ("\n\n");
00322     }
00323 
00324     /*printf ("=== "); fflush (stdout); 
00325     aux -> print (); printf (" := "); fflush (stdout);
00326     aux -> Image () -> print (); printf ("\n");*/
00327   }
00328 
00329   for (unsigned int i = iters2erase.size (); i--;)
00330     constraints_. erase (iters2erase [i]);
00331 
00332   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) { 
00333     // Use with caution. Bounds on auxs are not defined yet, so valgrind complains
00334     printf ("done with standardization: (careful, bounds below can be nonsense)\n");
00335     print (); 
00336   }
00337 
00338   delete auxSet_;
00339 
00340   // Create evaluation order ////////////////////////////////////////////////////
00341 
00342   // reallocate space for enlarged set of variables
00343   domain_.current () -> resize (nVars ());
00344 
00345   //graph_ -> print ();
00346   graph_ -> createOrder ();
00347   //graph_ -> print ();
00348 
00349   assert (graph_ -> checkCycles () == false);
00350 
00351   // Fill numbering structure /////////////////////////////////////////////////
00352 
00353   int n = nVars ();
00354   numbering_ = new int [n];
00355   std::set <DepNode *, compNode> vertices = graph_ -> Vertices ();
00356 
00357   for (std::set <DepNode *, compNode>::iterator i = vertices.begin ();
00358        i != vertices.end (); ++i)
00359 
00360     numbering_ [(*i) -> Order ()] = (*i) -> Index (); 
00361 
00363 
00364   // do initial bound propagation
00365 
00366   for (int i = 0; i < nVars (); i++) {
00367 
00368     int ord = numbering_ [i];
00369 
00370     if (variables_ [ord] -> Type () == AUX) {
00371 
00372       // initial auxiliary bounds are infinite (they are later changed
00373       // through branching)
00374 
00375       if (variables_ [ord] -> Index () >= nOrigVars_) { // and one that was not an original, originally...
00376 
00377         domain_.lb (ord) = -COIN_DBL_MAX;
00378         domain_.ub (ord) =  COIN_DBL_MAX;
00379       }
00380 
00381       //printf ("from "); variables_ [ord] -> Lb    () -> print (); 
00382 
00383       // tighten them with propagated bounds
00384       variables_ [ord] -> crossBounds ();
00385 
00386       //printf ("to "); variables_ [ord] -> Lb    () -> print (); printf (", now eval\n");
00387 
00388       if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00389         printf (":::: %3d %10g [%10g, %10g]", 
00390                 ord, domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
00391       }
00392 
00393       // and evaluate them
00394       domain_.x  (ord) = (*(variables_ [ord] -> Image ())) ();
00395       domain_.lb (ord) = (*(variables_ [ord] -> Lb    ())) ();
00396       domain_.ub (ord) = (*(variables_ [ord] -> Ub    ())) ();
00397 
00398       // if (DefVarDiffSet.find (ord) != DefVarDiffSet.end ()) {
00399 
00400       //        domain_.lb (ord) =
00401       //        domain_.ub (ord) = 0.;
00402       // }
00403 
00404       if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
00405         printf (" --> %10g [%10g, %10g] [", 
00406                 domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
00407         variables_ [ord] -> Lb () -> print (); printf (",");
00408         variables_ [ord] -> Ub () -> print (); printf ("]\n");
00409       }
00410 
00411       bool integer = variables_ [ord] -> isInteger ();
00412 
00413       if (integer) {
00414         domain_.lb (ord) = ceil  (domain_.lb (ord) - COUENNE_EPS);
00415         domain_.ub (ord) = floor (domain_.ub (ord) + COUENNE_EPS);
00416       }
00417     }
00418   }
00419 
00420   // TODO: resolve duplicate index in exprQuad before restoring this
00421 
00422   // remove duplicates
00423 
00424   std::string delete_redund;
00425 
00426   if (bonBase_) bonBase_ -> options () -> GetStringValue ("delete_redundant", delete_redund, "couenne."); 
00427   else  delete_redund = "yes";
00428 
00429   if (delete_redund == "yes")
00430 
00431     // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
00432     for (std::vector <exprVar *>::iterator i = variables_.begin (); 
00433        i != variables_.end (); ++i)
00434 
00435     if (((*i) -> Type () == AUX) && ((*i) -> sign () == expression::AUX_EQ)) {
00436 
00437       int type = (*i) -> Image () -> Type ();
00438 
00439       if ((type == VAR) || (type == AUX)) {
00440 
00441         // found w_k = x_h. 
00442         // 
00443         // Check if either is integer, the survivor will be integer too
00444         // Replace all occurrences of w_k with x_h
00445 
00446         /*printf ("redundancy: "); 
00447         (*i)             -> print (); printf (" := "); 
00448         (*i) -> Image () -> print (); printf ("\n");*/
00449 
00450         // use the info on the variable to be discarded: tighter
00451         // bounds and integrality that the replacement might not have.
00452 
00453         int 
00454           indStays  = (*i) -> Image () -> Index (), // index h
00455           indLeaves = (*i)             -> Index (); // index k
00456 
00457         if (indStays == indLeaves)  // very strange case, w_h = x_h
00458           continue;
00459 
00460         // do not swap! x_h could be in turn an auxiliary...
00461         //
00462         //if (indStays > indLeaves) 
00463         //{int swap = indStays; indStays = indLeaves; indLeaves = swap;} // swap
00464 
00465         exprVar 
00466           *varStays  = variables_ [indStays],
00467           *varLeaves = variables_ [indLeaves];
00468 
00469         // intersect features of the two variables (integrality, bounds)
00470 
00471         varStays -> lb () = varLeaves -> lb () = CoinMax (varStays -> lb (), varLeaves -> lb ());
00472         varStays -> ub () = varLeaves -> ub () = CoinMin (varStays -> ub (), varLeaves -> ub ());
00473 
00474         if (varStays  -> isInteger () ||
00475             varLeaves -> isInteger ()) {
00476 
00477           varStays -> lb () = ceil  (varStays -> lb ());
00478           varStays -> ub () = floor (varStays -> ub ());
00479 
00480           if (varStays -> Type () == AUX)
00481             varStays -> setInteger (true);
00482           else {
00483             //expression *old = varStays; // !!! leak
00484             variables_ [indStays] = varStays = new exprIVar (indStays, &domain_);
00485             auxiliarize (varStays); // replace it everywhere in the problem
00486             //delete old;
00487           }
00488         }
00489 
00490         auxiliarize (varLeaves, varStays); // now replace occurrences of w_k with x_h
00491 
00492         //if (varLeaves -> Index () >= nOrigVars_) // why check? It's not there anymore.
00493         varLeaves -> zeroMult (); // disable this variable
00494       }
00495     }
00496 
00502   for (int ii=0; ii < nVars (); ii++) {
00503 
00504     int i = numbering_ [ii];
00505 
00506     if ((Var (i) -> Multiplicity () > 0) &&
00507         (Var (i) -> Type () == AUX) &&
00508         (Var (i) -> Image () -> isInteger ()) &&
00509         (Var (i) -> sign () == expression::AUX_EQ))
00510       Var (i) -> setInteger (true);
00511 
00512     //if (Var (i) -> Multiplicity () == 0)
00513     //Lb (i) = Ub (i) = 0.;
00514   }
00515 
00516   // check how many multiplications there are 
00517 
00518 //   int nmul = 0;
00519 //   // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
00520 //   for (std::vector <exprVar *>::iterator i = variables_.begin (); 
00521 //        i != variables_.end (); ++i) {
00522 
00523 //     if ((*i) -> Type () != AUX ||
00524 //      (*i) -> Multiplicity () <= 0) 
00525 //       continue;
00526 
00527 //     expression *img = (*i) -> Image ();
00528 
00529 //     if (img -> code () != COU_EXPRMUL) 
00530 //       continue;
00531 
00532 //     expression **args = img -> ArgList ();
00533 
00534 //     if ((args [0] -> Type () == AUX ||
00535 //       args [0] -> Type () == VAR) &&
00536 //      (args [1] -> Type () == AUX ||
00537 //       args [1] -> Type () == VAR))
00538 //       nmul++;
00539 //   }
00540 
00541 //   printf ("\nMULS: %d/%d\n", nmul, variables_.size());
00542 
00543 //   exit (-1);
00544 
00545   // TODO: re-compute ranks
00546 
00547   delete [] commuted_;  commuted_ = NULL;
00548   delete    graph_;     graph_    = NULL;
00549 
00550   return retval;
00551 }

Generated on Thu Sep 22 03:06:00 2011 by  doxygen 1.4.7