/home/coin/SVN-release/OS-2.4.2/Couenne/src/problem/CouenneLPtightenBoundsCLP.cpp

Go to the documentation of this file.
00001 /* $Id: CouenneLPtightenBoundsCLP.cpp 490 2011-01-14 16:07:12Z pbelotti $
00002  *
00003  * Name:    CouenneLPtightenBoundsCLP.cpp
00004  * Authors: Pietro Belotti, Carnegie Mellon University
00005  * Purpose: adaptation of OsiClpSolverInterface::tightenBounds()
00006  *
00007  * (C) Carnegie-Mellon University, 2009.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CouennePrecisions.hpp"
00012 #include "CouenneProblem.hpp"
00013 #include "CouenneCutGenerator.hpp"
00014 #include "CouenneExprVar.hpp"
00015 
00016 namespace Couenne {
00017 
00018 //#define COIN_DEVELOP 4
00019 
00020 // Tighten bounds. Returns -1 if infeasible, otherwise number of
00021 // variables tightened.
00022 template <class T> 
00023 int CouenneSolverInterface<T>::tightenBoundsCLP (int lightweight) {
00024 
00025   // Copied from OsiClpSolverInterface::tightenBounds
00026 
00027   int
00028     numberRows    = T::getNumRows(),
00029     numberColumns = T::getNumCols(),
00030     iRow, iColumn;
00031 
00032   const double * columnUpper = T::getColUpper();
00033   const double * columnLower = T::getColLower();
00034   const double * rowUpper = T::getRowUpper();
00035   const double * rowLower = T::getRowLower();
00036 
00037   // Column copy of matrix
00038   const double * element = T::getMatrixByCol()->getElements();
00039   const int * row = T::getMatrixByCol()->getIndices();
00040   const CoinBigIndex * columnStart = T::getMatrixByCol()->getVectorStarts();
00041   const int * columnLength = T::getMatrixByCol()->getVectorLengths();
00042   const double *objective = T::getObjCoefficients() ;
00043 
00044   double direction = T::getObjSense();
00045   double * down = new double [numberRows];
00046 
00047   if (lightweight)
00048     return tightenBoundsCLP_Light (lightweight);
00049 
00050   // NOT LIGHTWEIGHT /////////////////////////////////////////////////////
00051 
00052   double * up = new double [numberRows];
00053   double * sum = new double [numberRows];
00054   int * type = new int [numberRows];
00055   CoinZeroN(down,numberRows);
00056   CoinZeroN(up,numberRows);
00057   CoinZeroN(sum,numberRows);
00058   CoinZeroN(type,numberRows);
00059   double infinity = T::getInfinity();
00060 
00061   for (iColumn=0;iColumn<numberColumns;iColumn++) {
00062     CoinBigIndex start = columnStart[iColumn];
00063     CoinBigIndex end = start + columnLength[iColumn];
00064     double lower = columnLower[iColumn];
00065     double upper = columnUpper[iColumn];
00066     if (lower==upper) {
00067       for (CoinBigIndex j=start;j<end;j++) {
00068         int iRow = row[j];
00069         double value = element[j];
00070         sum[iRow]+=2.0*fabs(value*lower);
00071         if ((type[iRow]&1)==0)
00072           down[iRow] += value*lower;
00073         if ((type[iRow]&2)==0)
00074           up[iRow] += value*lower;
00075       }
00076     } else {
00077       for (CoinBigIndex j=start;j<end;j++) {
00078         int iRow = row[j];
00079         double value = element[j];
00080         if (value>0.0) {
00081           if ((type[iRow]&1)==0) {
00082             if (lower!=-infinity) {
00083               down[iRow] += value*lower;
00084               sum[iRow]+=fabs(value*lower);
00085             } else {
00086               type[iRow] |= 1;
00087             }
00088           }
00089           if ((type[iRow]&2)==0) {
00090             if (upper!=infinity) {
00091               up[iRow] += value*upper;
00092               sum[iRow]+=fabs(value*upper);
00093             } else {
00094               type[iRow] |= 2;
00095             }
00096           }
00097         } else {
00098           if ((type[iRow]&1)==0) {
00099             if (upper!=infinity) {
00100               down[iRow] += value*upper;
00101               sum[iRow]+=fabs(value*upper);
00102             } else {
00103               type[iRow] |= 1;
00104             }
00105           }
00106           if ((type[iRow]&2)==0) {
00107             if (lower!=-infinity) {
00108               up[iRow] += value*lower;
00109               sum[iRow]+=fabs(value*lower);
00110             } else {
00111               type[iRow] |= 2;
00112             }
00113           }
00114         }
00115       }
00116     }
00117   }
00118 
00119   int nTightened = 0;
00120   double tolerance = 1.0e-6;
00121 
00122   for (iRow=0;iRow<numberRows;iRow++) {
00123     if ((type[iRow]&1)!=0)
00124       down[iRow]=-infinity;
00125     if (down[iRow]>rowUpper[iRow]) {
00126       if (down[iRow]>rowUpper[iRow]+tolerance+1.0e-8*sum[iRow]) {
00127         // infeasible
00128 #ifdef COIN_DEVELOP
00129         printf("infeasible on row %d\n",iRow);
00130 #endif
00131         nTightened=-1;
00132         break;
00133       } else {
00134         down[iRow]=rowUpper[iRow];
00135       }
00136     }
00137     if ((type[iRow]&2)!=0)
00138       up[iRow]=infinity;
00139     if (up[iRow]<rowLower[iRow]) {
00140       if (up[iRow]<rowLower[iRow]-tolerance-1.0e-8*sum[iRow]) {
00141         // infeasible
00142 #ifdef COIN_DEVELOP
00143         printf("infeasible on row %d\n",iRow);
00144 #endif
00145         nTightened=-1;
00146         break;
00147       } else {
00148         up[iRow]=rowLower[iRow];
00149       }
00150     }
00151   }
00152 
00153   if (nTightened)
00154     numberColumns = 0; // so will skip
00155 
00156   for (iColumn=0;iColumn<numberColumns;iColumn++) {
00157     double lower = columnLower[iColumn];
00158     double upper = columnUpper[iColumn];
00159     double gap = upper-lower;
00160 
00161     if (!gap)
00162       continue;
00163 
00164     int canGo=0;
00165 
00166     CoinBigIndex
00167       start =         columnStart  [iColumn],
00168       end   = start + columnLength [iColumn];
00169 
00170     if (lower < -1.0e8 && upper > 1.0e8)
00171       continue; // Could do severe damage to accuracy
00172 
00173 
00174     // there was an ifInteger condition here. We do like tightened
00175     // bounds for continuous variables too, so we don't test for
00176     // integrality.
00177 
00178     std::vector <exprVar *> &vars = cutgen_ -> Problem () -> Variables ();
00179 
00180     {
00181       if (vars [iColumn] -> isInteger ()) {
00182 
00183         if (lower < ceil (lower - COUENNE_EPS) - COUENNE_EPS) {
00184 #ifdef COIN_DEVELOP
00185           printf("increasing lower bound on %d from %e to %e\n",iColumn,
00186                  lower,ceil(lower - COUENNE_EPS));
00187 #endif
00188           lower=ceil(lower - COUENNE_EPS);
00189           gap=upper-lower;
00190           T::setColLower(iColumn,lower);
00191         }
00192 
00193         if (upper > floor(upper + COUENNE_EPS) + COUENNE_EPS) {
00194 #ifdef COIN_DEVELOP
00195           printf("decreasing upper bound on %d from %e to %e\n",iColumn,
00196                  upper,floor(upper + COUENNE_EPS));
00197 #endif
00198           upper=floor(upper + COUENNE_EPS);
00199           gap=upper-lower;
00200           T::setColUpper(iColumn,upper);
00201         }
00202       }
00203 
00204       double newLower=lower;
00205       double newUpper=upper;
00206 
00207       for (CoinBigIndex j=start;j<end;j++) {
00208         int iRow = row[j];
00209         double value = element[j];
00210         if (value>0.0) {
00211           if ((type[iRow]&1)==0) {
00212             // has to be at most something
00213             if (down[iRow] + value*gap > rowUpper[iRow]+tolerance) {
00214               double newGap = (rowUpper[iRow]-down[iRow])/value;
00215               // adjust
00216               newGap += 1.0e-10*sum[iRow];
00217               if (vars [iColumn] -> isInteger ())
00218                 newGap = floor(newGap);
00219               if (lower+newGap<newUpper)
00220                 newUpper=lower+newGap;
00221             }
00222           }
00223           if (down[iRow]<rowLower[iRow])
00224             canGo |=1; // can't go down without affecting result
00225           if ((type[iRow]&2)==0) {
00226             // has to be at least something
00227             if (up[iRow] - value*gap < rowLower[iRow]-tolerance) {
00228               double newGap = (up[iRow]-rowLower[iRow])/value;
00229               // adjust
00230               newGap += 1.0e-10*sum[iRow];
00231               if (vars [iColumn] -> isInteger ())
00232                 newGap = floor(newGap);
00233               if (upper-newGap>newLower)
00234                 newLower=upper-newGap;
00235             }
00236           }
00237           if (up[iRow]>rowUpper[iRow])
00238             canGo |=2; // can't go up without affecting result
00239         } else {
00240           if ((type[iRow]&1)==0) {
00241             // has to be at least something
00242             if (down[iRow] - value*gap > rowUpper[iRow]+tolerance) {
00243               double newGap = -(rowUpper[iRow]-down[iRow])/value;
00244               // adjust
00245               newGap += 1.0e-10*sum[iRow];
00246               if (vars [iColumn] -> isInteger ())
00247                 newGap = floor(newGap);
00248               if (upper-newGap>newLower)
00249                 newLower=upper-newGap;
00250             }
00251           }
00252           if (up[iRow]>rowUpper[iRow])
00253             canGo |=1; // can't go down without affecting result
00254           if ((type[iRow]&2)==0) {
00255             // has to be at most something
00256             if (up[iRow] + value*gap < rowLower[iRow]-tolerance) {
00257               double newGap = -(up[iRow]-rowLower[iRow])/value;
00258               // adjust
00259               newGap += 1.0e-10*sum[iRow];
00260               if (vars [iColumn] -> isInteger ())
00261                 newGap = floor(newGap);
00262               if (lower+newGap<newUpper)
00263                 newUpper=lower+newGap;
00264             }
00265           }
00266           if (down[iRow]<rowLower[iRow])
00267             canGo |=2; // can't go up without affecting result
00268         }
00269       }
00270 
00271       if (newUpper<upper || newLower>lower) {
00272         nTightened++;
00273         if (newLower>newUpper) {
00274           // infeasible
00275 #if COIN_DEVELOP>1
00276           printf("infeasible on column %d\n",iColumn);
00277 #endif
00278           nTightened=-1;
00279           break;
00280         } else {
00281           T::setColLower(iColumn,newLower);
00282           T::setColUpper(iColumn,newUpper);
00283         }
00284         for (CoinBigIndex j=start;j<end;j++) {
00285           int iRow = row[j];
00286           double value = element[j];
00287           if (value>0.0) {
00288             if ((type[iRow]&1)==0) down [iRow] += value*(newLower-lower);
00289             if ((type[iRow]&2)==0) up   [iRow] += value*(newUpper-upper);
00290           } else {
00291             if ((type[iRow]&1)==0) down [iRow] += value*(newUpper-upper);
00292             if ((type[iRow]&2)==0) up   [iRow] += value*(newLower-lower);
00293           }
00294         }
00295       } else {
00296 
00297         if (canGo!=3) {
00298 
00299           double objValue = direction*objective[iColumn];
00300 
00301           if (objValue>=0.0&&(canGo&1)==0) {
00302 #if COIN_DEVELOP>2
00303             printf("dual fix down on column %d\n",iColumn);
00304 #endif
00305             nTightened++;
00306             T::setColUpper(iColumn,lower);
00307           } else if (objValue<=0.0 && (canGo&2)==0) {
00308 #if COIN_DEVELOP>2
00309             printf("dual fix up on column %d\n",iColumn);
00310 #endif
00311             nTightened++;
00312             T::setColLower(iColumn,upper);
00313           }
00314         }           
00315       }
00316     }
00317 
00318 //     else {
00319 
00320 //       // CONTINUOUS //////////////////////////////////////////
00321 
00322 //       // just do dual tests
00323 //       for (CoinBigIndex j=start;j<end;j++) {
00324 //      int iRow = row[j];
00325 //      double value = element[j];
00326 //      if (value>0.0) {
00327 //        if (down [iRow] < rowLower [iRow]) canGo |=1; // can't go down without affecting result
00328 //        if (up   [iRow] > rowUpper [iRow]) canGo |=2; // can't go up   without affecting result
00329 //      } else {
00330 //        if (up   [iRow] > rowUpper [iRow]) canGo |=1; // can't go down without affecting result
00331 //        if (down [iRow] < rowLower [iRow]) canGo |=2; // can't go up   without affecting result
00332 //      }
00333 //       }
00334 
00335 //       if (canGo!=3) {
00336 //      double objValue = direction*objective[iColumn];
00337 //      if (objValue>=0.0&&(canGo&1)==0) {
00338 // #if COIN_DEVELOP>2
00339 //        printf("dual fix down on continuous column %d lower %g\n",
00340 //               iColumn,lower);
00341 // #endif
00342 //        // Only if won't cause numerical problems
00343 //        if (lower>-1.0e10) {
00344 //          nTightened++;;
00345 //          setColUpper(iColumn,lower);
00346 //        }
00347 //      } else if (objValue<=0.0&&(canGo&2)==0) {
00348 // #if COIN_DEVELOP>2
00349 //        printf("dual fix up on continuous column %d upper %g\n",
00350 //               iColumn,upper);
00351 // #endif
00352 //        // Only if won't cause numerical problems
00353 //        if (upper<1.0e10) {
00354 //          nTightened++;;
00355 //          setColLower(iColumn,upper);
00356 //        }
00357 //      }
00358 //       }
00359 //     }
00360 
00361   }
00362 
00363   delete [] type;
00364   delete [] down;
00365   delete [] up;
00366   delete [] sum;
00367 
00368   return nTightened;
00369 }
00370   
00371 }

Generated on Wed Nov 30 03:04:07 2011 by  doxygen 1.4.7