/home/coin/SVN-release/OS-2.1.0/Couenne/src/disjunctive/getDisjunctions.cpp

Go to the documentation of this file.
00001 /* $Id: getDisjunctions.cpp 236 2009-07-19 16:15:55Z pbelotti $ */
00002 /*
00003  * Name:    getDisjunctions.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: get nonlinear (and integer?) disjunctions for the problem
00006  *
00007  * (C) Carnegie-Mellon University, 2008. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CouenneCutGenerator.hpp"
00012 #include "CouenneProblem.hpp"
00013 #include "CouenneObject.hpp"
00014 #include "CouenneBranchingObject.hpp"
00015 #include "CouenneDisjCuts.hpp"
00016 
00017 
00019 int CouenneDisjCuts::getDisjunctions (std::vector <std::pair <OsiCuts *, OsiCuts *> > &disjs, 
00020                                       OsiSolverInterface &si, 
00021                                       OsiCuts &cs, 
00022                                       const CglTreeInfo &info) const {
00023 
00024   OsiBranchingInformation brInfo (&si, 
00025                                   true,   // bool normalSolver,
00026                                   false); // bool copySolution=false);
00027 
00028   if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS))
00029     for (int i=0; i<si.getNumCols (); i++)
00030       printf ("%3d. %8g  [%8g %8g]\n", i, 
00031               brInfo. solution_ [i], 
00032               brInfo. lower_    [i], 
00033               brInfo. upper_    [i]);
00034 
00035   // get list of best disjunction as though we were branching
00036   branchingMethod_ -> setupList (&brInfo, true); // initialize
00037 
00038   int 
00039     ncols  = si.getNumCols (),
00040     retval = COUENNE_FEASIBLE;
00041 
00042   int nobjs = branchingMethod_ -> numberUnsatisfied ();
00043 
00044   if (nobjs) {
00045 
00046     if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) 
00047       printf ("---   %d disjunctive objects\n", nobjs);
00048 
00049     const int *candidates = branchingMethod_ -> candidates ();
00050     OsiObject **objs = si. objects ();
00051 
00052     // enumerate first (most infeasible, or most promising) objects of
00053     // the list
00054 
00055     for (int candInd = 0,      num = 0; 
00056          (candInd < nobjs) && (num < numDisjunctions_); 
00057          candInd++) {
00058 
00059       OsiObject *Object = objs [candidates [candInd]];
00060       CouenneObject *cObj = dynamic_cast <CouenneObject *> (Object);
00061 
00062       if (cObj &&                                          // IF    a nonlinear object BUT:
00063           ((cObj -> checkInfeasibility (&brInfo) == 0.) || //       not violated
00064            (cObj -> isCuttable ())))                       //    or we are on the "good" (convex) side
00065         continue;                                          // THEN skip
00066 
00067       bool     feasLeft = true,  feasRight = true;
00068       OsiCuts *leftCuts = NULL, *rightCuts = NULL;
00069 
00070       const double 
00071         *saveLower = CoinCopyOfArray (si.getColLower (), ncols),
00072         *saveUpper = CoinCopyOfArray (si.getColUpper (), ncols);
00073 
00074       OsiBranchingObject *brObj = Object -> createBranch (&si, &brInfo, 0); // down!
00075 
00076       if (jnlst_ -> ProduceOutput (J_VECTOR, J_DISJCUTS))
00077         printf ("---   cand [%d] is %d: x_%d <>= %g [%g,%g]\n", 
00078               candInd, 
00079                 candidates [candInd],
00080                 dynamic_cast <CouenneBranchingObject *> (brObj) -> variable () -> Index (), 
00081                 brObj -> value (),
00082                 couenneCG_->Problem()->Lb(dynamic_cast<CouenneBranchingObject*>(brObj)
00083                                           ->variable()->Index ()),
00084                 couenneCG_->Problem()->Ub(dynamic_cast<CouenneBranchingObject*>(brObj)
00085                                           ->variable()->Index ()));
00086 
00087       // left branch ////////////////////////////////////////////////////////////
00088 
00089       if ((brObj -> branch (&si) >= COUENNE_INFINITY) ||
00090           // Bound tightening if not a CouenneObject -- explicit 
00091           (!cObj && !BranchingFBBT (couenneCG_ -> Problem (), Object, &si)))
00092         feasLeft = false;                        // left subtree infeasible
00093       else leftCuts = getSingleDisjunction (si); // feasible, store new bounds in OsiCuts
00094 
00095       // restore initial bounds
00096       si.setColLower (saveLower);
00097       si.setColUpper (saveUpper);
00098 
00099       delete brObj;
00100 
00101       // right branch ////////////////////////////////////////////////////////////
00102 
00103       brObj = Object -> createBranch (&si, &brInfo, 1); // up!
00104 
00105       if ((brObj -> branch (&si) >= COUENNE_INFINITY)  ||
00106           // Bound tightening if not a CouenneObject -- explicit 
00107           (!cObj && !BranchingFBBT (couenneCG_ -> Problem (), Object, &si)))
00108         feasRight = false;                        // right subtree infeasible
00109       else rightCuts = getSingleDisjunction (si); // feasible, store new bounds in OsiCuts
00110 
00111       delete brObj;
00112 
00113       // restore initial bounds
00114       si.setColLower (saveLower);
00115       si.setColUpper (saveUpper);
00116 
00117       // done with generating a disjunction. Is any (or both) side infeasible?
00118 
00119       if (feasLeft && feasRight) {
00120 
00121         // feasible, which variables tightened?
00122 
00123         std::pair <OsiCuts *, OsiCuts *> newpair;
00124         newpair.first  = leftCuts;
00125         newpair.second = rightCuts;
00126 
00127         disjs.push_back (newpair);
00128 
00129         num++;
00130 
00131       } else if (feasLeft) {
00132 
00133         jnlst_ -> Printf (J_VECTOR, J_DISJCUTS, "---   disj: infeasible right\n");
00134         applyColCuts (si, leftCuts); // !!! culprit
00135         retval = COUENNE_TIGHTENED;
00136 
00137       } else if (feasRight) {
00138 
00139         jnlst_ -> Printf (J_VECTOR, J_DISJCUTS, "---   infeasible left\n");
00140         applyColCuts (si, rightCuts); // !!! culprit
00141         retval = COUENNE_TIGHTENED;
00142 
00143       } else {
00144         jnlst_ -> Printf (J_VECTOR, J_DISJCUTS, "---   infeasible!!!\n");
00145         retval = COUENNE_INFEASIBLE;
00146         break;
00147       }
00148 
00149       delete [] saveLower;
00150       delete [] saveUpper;
00151     }
00152   }
00153 
00154   // disjunctions (of bounding boxes) generated. Do some extra
00155   // tightening that may result from tightened sides of the
00156   // disjunctions
00157 
00158   if (retval != COUENNE_INFEASIBLE) {
00159 
00160     if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS))
00161       printf ("have %d disjunctions\n", disjs.size ());
00162 
00163     // sanity check: for each disjunction, check if left and right
00164     // part are feasible with (possibly improved) bounds of si --
00165     // these may have tightened and now exclude either (or both!)
00166     // sides of a disjunction
00167 
00168     for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjs.begin ();
00169          disjI != disjs.end (); ++disjI) {
00170 
00171       if   (checkDisjSide (si, disjI -> first)  == COUENNE_INFEASIBLE) { // left  side infeasible?
00172         if (checkDisjSide (si, disjI -> second) == COUENNE_INFEASIBLE) { // right side infeasible?
00173 
00174           retval = COUENNE_INFEASIBLE;
00175           break;
00176 
00177         } else {
00178 
00179           if (jnlst_ -> ProduceOutput (J_VECTOR, J_DISJCUTS))
00180             printf ("---   infeasible left [2]!\n");
00181           applyColCuts (si, disjI -> second);
00182           retval = COUENNE_TIGHTENED;
00183         }
00184       } else if (checkDisjSide (si, disjI -> second) == COUENNE_INFEASIBLE) { // right side
00185 
00186         if (jnlst_ -> ProduceOutput (J_VECTOR, J_DISJCUTS))
00187           printf ("---   infeasible right [2]!\n");
00188         applyColCuts (si, disjI -> first);
00189         retval = COUENNE_TIGHTENED;
00190       }
00191     }
00192 
00193     if (retval == COUENNE_INFEASIBLE)
00194       return retval;
00195 
00196     // merge tightened disjunction by intersecting si's bounding box
00197     // with intersections of smallest boxes, one per disjunction,
00198     // containing each both sides of the disjunction
00199 
00200     for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjs.begin ();
00201          disjI != disjs.end (); ++disjI) {
00202 
00203       CoinPackedVector lowerChg, upperChg;
00204 
00205       // find smallest box containing two disjunctions
00206       getBoxUnion (si, disjI -> first, disjI -> second, lowerChg, upperChg);
00207 
00208       if ((lowerChg.getNumElements () > 0) || 
00209           (upperChg.getNumElements () > 0)) {
00210 
00211         OsiColCut cut;
00212         cut.setLbs (lowerChg);
00213         cut.setUbs (upperChg);
00214         applyColCuts (si, &cut);
00215 
00216         cs.insert (cut);
00217       }
00218     }   
00219   }
00220 
00221   return retval;
00222 }
00223 
00224 
00226 void CouenneDisjCuts::applyColCuts (OsiSolverInterface &si, OsiCuts *cuts) const {
00227 
00228   if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS)) {
00229     printf ("applying cuts to SI:\n");
00230     // apply column cuts to si
00231     for (int i = cuts -> sizeColCuts (); i--;)
00232       cuts -> colCutPtr (i) -> print ();
00233     printf ("--------------------\n");
00234   }
00235 
00236   // apply column cuts to si
00237   for (int i = cuts -> sizeColCuts (); i--;)
00238     applyColCuts (si, cuts -> colCutPtr (i));
00239 }
00240 
00241 
00243 void CouenneDisjCuts::applyColCuts (OsiSolverInterface &si, OsiColCut *cut) const {
00244 
00245   // apply single column cut to si
00246   if (jnlst_ -> ProduceOutput (J_VECTOR, J_DISJCUTS)) {
00247     printf ("tightening bounds: "); 
00248     cut -> print ();
00249   }
00250 
00251   const CoinPackedVector
00252     &lbs = cut -> lbs (),
00253     &ubs = cut -> ubs ();
00254 
00255   const int    *lind = lbs.getIndices  (), *uind = ubs.getIndices  ();
00256   const double *lval = lbs.getElements (), *uval = ubs.getElements (),
00257     *oldLower = si.getColLower (), 
00258     *oldUpper = si.getColUpper ();
00259 
00260   for (int j=lbs.getNumElements (); j--; lval++, lind++)
00261     if (*lval > oldLower [*lind] + COUENNE_EPS) 
00262       si.setColLower (*lind, *lval);
00263 
00264   for (int j=ubs.getNumElements (); j--; uval++, uind++)
00265     if (*uval < oldUpper [*uind] - COUENNE_EPS) 
00266       si.setColUpper (*uind, *uval);
00267 }

Generated on Tue Mar 30 03:04:37 2010 by  doxygen 1.4.7