/home/coin/SVN-release/OS-2.4.1/Couenne/src/bound_tightening/FixPointGenCuts.cpp

Go to the documentation of this file.
00001 /* $Id: FixPointGenCuts.cpp 694 2011-06-18 20:13:17Z stefan $
00002  *
00003  * Name:    FixPointGenCuts.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: Fix point bound tightener -- separator
00006  *
00007  * (C) Pietro Belotti, 2010.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CoinHelperFunctions.hpp"
00012 #include "OsiClpSolverInterface.hpp"
00013 #include "OsiCuts.hpp"
00014 #include "CoinTime.hpp"
00015 
00016 #include "CouenneFixPoint.hpp"
00017 
00018 #include "CouenneProblem.hpp"
00019 #include "CouennePrecisions.hpp"
00020 #include "CouenneExprVar.hpp"
00021 #include "CouenneInfeasCut.hpp"
00022 
00023 using namespace Ipopt;
00024 using namespace Couenne;
00025 
00027 
00028 void CouenneFixPoint::generateCuts (const OsiSolverInterface &si,
00029                                     OsiCuts &cs,
00030                                     const CglTreeInfo treeInfo) const {
00031 
00040 
00041   if (firstCall_) 
00042     firstCall_ = false;
00043   else
00044     if (!(problem_ -> fbbtReachedIterLimit ()))
00045       return;
00046 
00047   if (isWiped (cs))
00048     return;
00049 
00050   if (treeInfo.inTree && 
00051       treeInfo.level > 0 &&
00052       treeInfo.pass > 1)
00053     return;
00054 
00055   int nInitTightened = nTightened_;
00056 
00057   if (treeInfo.inTree && 
00058       treeInfo.level <= 0) {
00059 
00060     problem_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "Fixed Point FBBT: "); 
00061     fflush (stdout);
00062   }
00063 
00064   ++nRuns_;
00065 
00066   double now = CoinCpuTime ();
00067 
00068   problem_ -> domain () -> push (&si, &cs);
00069 
00088 
00089   OsiSolverInterface *fplp = NULL;
00090 
00091   if (true) { // placeholder for later selection of LP solver among
00092               // those available
00093 
00094     fplp = new OsiClpSolverInterface;
00095   }
00096 
00113 
00115 
00116   const CoinPackedMatrix *A = si. getMatrixByRow ();
00117 
00118   const int
00119        n     = si.  getNumCols  (),
00120        m     = si.  getNumRows  (),
00121        nCuts = cs.sizeRowCuts   (),
00122     *ind     = A -> getIndices  ();
00123 
00124   const double
00125     *lb  = si.  getColLower (),
00126     *ub  = si.  getColUpper (),
00127     *rlb = si.  getRowLower (),
00128     *rub = si.  getRowUpper (),
00129     *coe = A -> getElements ();
00130 
00131   if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING))
00132     for (int i=0; i<n; i++) 
00133       printf ("----------- x_%d in [%g,%g]\n", i, lb [i], ub [i]);
00134 
00135   // turn off logging
00136   fplp -> messageHandler () -> setLogLevel (0);
00137 
00138   // add lvars and uvars to the new problem
00139   for (int i=0; i<n; i++) {
00140     bool isActive = problem_ -> Var (i) -> Multiplicity () > 0;
00141     fplp -> addCol (0, NULL, NULL, lb [i], ub [i], isActive ? -1. : 0.); // xL_i
00142   }
00143 
00144   for (int i=0; i<n; i++) {
00145     bool isActive = problem_ -> Var (i) -> Multiplicity () > 0;
00146     fplp -> addCol (0, NULL, NULL, lb [i], ub [i], isActive ? +1. : 0.); // xU_i
00147   }
00148 
00149   if (extendedModel_) {
00150 
00151     for (int j=0; j<m; j++) fplp -> addCol (0, NULL, NULL, rlb [j],       COIN_DBL_MAX, 0.); // bL_j
00152     for (int j=0; j<m; j++) fplp -> addCol (0, NULL, NULL, -COIN_DBL_MAX, rub [j],      0.); // bU_j
00153   }
00154 
00155   // Scan each row of the matrix 
00156   
00157   for (int j=0; j<m; j++) { // for each row
00158 
00159     int nEl = A -> getVectorSize (j); // # elements in each row
00160 
00161     if (!nEl)
00162       continue;
00163 
00164     if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING)) {
00165 
00166       printf ("row %4d, %4d elements: ", j, nEl);
00167 
00168       for (int i=0; i<nEl; i++) {
00169         printf ("%+g x%d ", coe [i], ind [i]);
00170         fflush (stdout);
00171       }
00172 
00173       printf ("\n");
00174     }
00175 
00176     // create cuts for the xL and xU elements //////////////////////
00177 
00178     if (extendedModel_ || rlb [j] > -COUENNE_INFINITY) 
00179       for (int i=0; i<nEl; i++) 
00180         createRow (-1, ind [i], n, fplp, ind, coe, rlb [j], nEl, extendedModel_, j, m + nCuts); // downward constraints -- on x_i
00181 
00182     if (extendedModel_ || rub [j] <  COUENNE_INFINITY) 
00183       for (int i=0; i<nEl; i++) 
00184         createRow (+1, ind [i], n, fplp, ind, coe, rub [j], nEl, extendedModel_, j, m + nCuts); // downward constraints -- on x_i
00185 
00186     // create (at most 2) cuts for the bL and bU elements //////////////////////
00187 
00188     if (extendedModel_) {
00189       createRow (-1, 2*n     + j, n, fplp, ind, coe, rlb [j], nEl, extendedModel_, j, m + nCuts); // upward constraints -- on bL_i
00190       createRow (+1, 2*n + m + j, n, fplp, ind, coe, rub [j], nEl, extendedModel_, j, m + nCuts); // upward constraints -- on bU_i
00191     }
00192 
00193     ind += nEl;
00194     coe += nEl;
00195   }
00196 
00197   // similarly, scan previous cuts in cs //////////////////////////////////////
00198 
00199   for (int j = 0, jj = nCuts; jj--; j++) {
00200 
00201     // create cuts for the xL and xU elements //////////////////////
00202 
00203     OsiRowCut *cut = cs.rowCutPtr (j);
00204 
00205     const CoinPackedVector &row = cut -> row ();
00206 
00207     const int 
00208       nEl = row.getNumElements (),
00209       *ind = row.getIndices ();
00210 
00211     const double *coe = row.getElements ();
00212 
00213     if (extendedModel_) {
00214       fplp -> addCol (0, NULL, NULL, cut -> lb (),       COIN_DBL_MAX, 0.); // bL_j
00215       fplp -> addCol (0, NULL, NULL, -COIN_DBL_MAX, cut -> ub (),      0.); // bU_j
00216     }
00217 
00218     if (extendedModel_ || cut -> lb () > -COUENNE_INFINITY) 
00219       for (int i=0; i<nEl; i++) 
00220         createRow (-1, ind [i], n, fplp, ind, coe, cut -> lb (), nEl, extendedModel_, m + j, m + nCuts); // downward constraints -- on x_i
00221 
00222     if (extendedModel_ || cut -> ub () <  COUENNE_INFINITY) 
00223       for (int i=0; i<nEl; i++) 
00224         createRow (+1, ind [i], n, fplp, ind, coe, cut -> ub (), nEl, extendedModel_, m + j, m + nCuts); // downward constraints -- on x_i
00225 
00226     // create (at most 2) cuts for the bL and bU elements
00227 
00228     if (extendedModel_) {
00229       createRow (-1, 2*n             + j, n, fplp, ind, coe, cut -> lb (), nEl, extendedModel_, m + j, m + nCuts); // upward constraints -- on bL_i
00230       createRow (+1, 2*n + m + nCuts + j, n, fplp, ind, coe, cut -> ub (), nEl, extendedModel_, m + j, m + nCuts); // upward constraints -- on bU_i
00231     }
00232 
00233     ind += nEl;
00234     coe += nEl;
00235   }
00236 
00237   // finally, add consistency cuts, bL <= bU
00238 
00239   if (extendedModel_)
00240 
00241     for (int j=0; j<m; j++) { // for each row
00242 
00243       int    ind [2] = {2*n + j, 2*n + m + j};
00244       double coe [2] = {1.,      -1.};
00245 
00246       CoinPackedVector row (2, ind, coe);
00247       fplp -> addRow (row, -COIN_DBL_MAX, 0.);
00248     }
00249 
00252 
00253   fplp -> setObjSense (-1.); // we want to maximize 
00254 
00255   //printf ("(writing lp) ");
00256   //fplp -> writeLp ("fplp");
00257 
00258   fplp -> initialSolve ();
00259 
00260   if (fplp -> isProvenOptimal ()) {
00261 
00262     // if problem not solved to optimality, bounds are useless
00263 
00264     const double 
00265       *newLB = fplp -> getColSolution (),
00266       *newUB = newLB + n,
00267       *oldLB = si. getColLower (),
00268       *oldUB = si. getColUpper ();
00269 
00270     // check old and new bounds
00271 
00272     int 
00273       *indLB = new int [n],
00274       *indUB = new int [n],
00275       ntightenedL = 0,
00276       ntightenedU = 0;
00277 
00278     double 
00279       *valLB = new double [n],
00280       *valUB = new double [n];
00281 
00282     for (int i=0; i<n; i++) {
00283 
00284       if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING))
00285         printf ("x%d: [%g,%g] --> [%g,%g]\n", i, 
00286                 oldLB [i], oldUB [i], 
00287                 newLB [i], newUB [i]);
00288 
00289       if (newLB [i] > oldLB [i] + COUENNE_EPS) {
00290 
00291         indLB [ntightenedL]   = i;
00292         valLB [ntightenedL++] = newLB [i];
00293 
00294         ++nTightened_;
00295       }
00296 
00297       if (newUB [i] < oldUB [i] - COUENNE_EPS) {
00298 
00299         indUB [ntightenedU]   = i;
00300         valUB [ntightenedU++] = newUB [i];
00301 
00302         ++nTightened_;
00303       }
00304     }
00305 
00306     if (ntightenedL || ntightenedU) {
00307 
00308       OsiColCut newBound;
00309 
00310       newBound.setLbs (ntightenedL, indLB, valLB);
00311       newBound.setUbs (ntightenedU, indUB, valUB);
00312 
00313       cs.insert (newBound);
00314     }
00315 
00316     delete [] indLB;
00317     delete [] indUB;
00318     delete [] valLB;
00319     delete [] valUB;
00320 
00321     CPUtime_ += CoinCpuTime () - now;
00322 
00323     if (treeInfo.inTree && 
00324         treeInfo.level <= 0) {
00325 
00326       problem_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "%d bounds tightened (%g seconds)\n", 
00327                                       nTightened_ - nInitTightened, CoinCpuTime () - now); 
00328     }
00329 
00330   } else
00331     if (treeInfo.inTree && 
00332         treeInfo.level <= 0)
00333       problem_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, " FPLP infeasible or unbounded.\n");
00334 
00335   delete fplp;
00336 
00337   problem_ -> domain () -> pop ();
00338 }
00339 
00340 
00341 // single cut creation. Parameters:
00342 //
00343 //  1) sign:     tells us whether this is a <= or a >= (part of a) constraint.
00344 //  2) indexVar: index of variable we want to do upward or downward bt
00345 //  3) nVars:    number of variables in the original problems (original +
00346 //               auxiliaries). Used to understand if we are adding an
00347 //               up or a down constraint
00348 //  4) p:        solver interface to which we are adding constraints
00349 //  5) indices:  vector containing indices of the linearization constraint (the    i's)
00350 //  6) coe:                        coeffs                                       a_ji's
00351 //  7) rhs:      right-hand side of constraint
00352 //  8) nEl:      number of elements of this linearization cut
00353 //  9) extMod:   extendedModel_
00354 // 10) indCon:   index of constraint being treated (and corresponding bL, bU)
00355 // 11) nCon:     number of constraints
00356 
00357 void CouenneFixPoint::createRow (int sign,
00358                                  int indexVar,
00359                                  int nVars,
00360                                  OsiSolverInterface *p,
00361                                  const int *indices,
00362                                  const double *coe,
00363                                  const double rhs,
00364                                  const int nEl,
00365                                  bool extMod,
00366                                  int indCon,
00367                                  int nCon) const {
00368 
00470 
00471 
00472   if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING)) {
00473 
00474     printf ("creating constraint from: ");
00475 
00476     for (int i=0; i<nEl; i++)
00477       printf ("%+g x%d ", coe [i], indices [i]);
00478 
00479     printf ("%c= %g for variable index %d: ", sign > 0 ? '<' : '>', rhs, indexVar);
00480   }
00481 
00482   int nTerms = nEl;
00483 
00484   if (extMod) 
00485     nTerms++; // always add one element when using extended model
00486 
00487   int    *iInd = new int    [nTerms];
00488   double *elem = new double [nTerms];
00489 
00490   // coefficients are really easy
00491 
00492   CoinCopyN (coe, nEl, elem);
00493 
00494   if (extMod) {
00495     elem [nEl] = -1.; // extended model, coefficient for bL or bU
00496     iInd [nEl] = 2*nVars + indCon + ((sign > 0) ? nCon : 0);
00497   }
00498 
00499   // indices are not so easy...
00500 
00501   for (int k=0; k<nEl; k++) {
00502 
00503     int curInd = indices [k];
00504 
00505     iInd [k] = curInd; // Begin with xL_i, same index as x_i in the
00506                        // original model. Should add n depending on a
00507                        // few things... 
00508 
00509     if (curInd == indexVar) { // x_k is x_i itself
00510       if (((sign > 0) && (coe [k] > 0.)) || 
00511           ((sign < 0) && (coe [k] < 0.)))
00512 
00513       iInd [k] += nVars;
00514 
00515     } else if (((coe [k] > 0.) && (sign < 0)) ||
00516                ((coe [k] < 0.) && (sign > 0)))
00517 
00518       iInd [k] += nVars;
00519   }
00520 
00521   CoinPackedVector vec (nTerms, iInd, elem);
00522 
00523   double 
00524     lb = sign > 0 ? -COIN_DBL_MAX : extMod ? 0. : rhs,
00525     ub = sign < 0 ? +COIN_DBL_MAX : extMod ? 0. : rhs;
00526 
00527   p -> addRow (vec, lb, ub);
00528 
00529   // Update time spent doing this
00530 
00531   if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_BOUNDTIGHTENING)) {
00532 
00533     for (int i=0; i<nEl; i++)
00534       printf ("%+g x%d ", elem [i], iInd [i]);
00535 
00536     printf ("in [%g,%g]\n", lb, ub);
00537   }
00538 
00539   // OsiRowCut *cut = new OsiRowCut (lb, ub, nTerms, nTerms, iInd, elem);
00540   // cut -> print ();
00541   // delete cut;
00542 
00543   delete [] iInd;
00544   delete [] elem;
00545 }

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