#include <CglSimpleRounding.hpp>
Inheritance diagram for CglSimpleRounding:
Public Member Functions | |
Generate Cuts | |
virtual void | generateCuts (const OsiSolverInterface &si, OsiCuts &cs) const |
Constructors and destructors | |
CglSimpleRounding () | |
Default constructor. | |
CglSimpleRounding (const CglSimpleRounding &) | |
Copy constructor. | |
CglSimpleRounding & | operator= (const CglSimpleRounding &rhs) |
Assignment operator. | |
virtual | ~CglSimpleRounding () |
Destructor. | |
Private Member Functions | |
Private methods | |
bool | deriveAnIntegerRow (const OsiSolverInterface &si, int rowIndex, const CoinShallowPackedVector &matrixRow, CoinPackedVector &irow, double &b, bool *negative) const |
Derive a <= inequality in integer variables from the rowIndex-th constraint. | |
int | power10ToMakeDoubleAnInt (int size, const double *x, double dataTol) const |
Greatest common denominators methods | |
int | gcd (int a, int b) const |
Returns the greatest common denominator of two positive integers, a and b. | |
int | gcdv (int n, const int *const vi) const |
Private Attributes | |
Private member data | |
double | epsilon_ |
A value within an epsilon_ neighborhood of 0 is considered to be 0. | |
Friends | |
void | CglSimpleRoundingUnitTest (const OsiSolverInterface *siP, const std::string mpdDir) |
This class generates simple rounding cuts via the following method: For each contraint, attempt to derive a <= inequality in all integer variables by netting out any continuous variables. Divide the resulting integer inequality through by the greatest common denomimator (gcd) of the lhs coefficients. Round down the rhs.
Warning: Use with careful attention to data precision.
(Reference: Nemhauser and Wolsey, Integer and Combinatorial Optimization, 1988, pg 211.)
Definition at line 26 of file CglSimpleRounding.hpp.
|
Returns the greatest common denominator of a vector of positive integers, vi, of length n. Definition at line 141 of file CglSimpleRounding.hpp. References gcd(). Referenced by generateCuts().
|
|
Generate simple rounding cuts for the model accessed through the solver interface. Insert generated cuts into the cut set cs. Implements CglCutGenerator. Definition at line 20 of file CglSimpleRounding.cpp. References deriveAnIntegerRow(), epsilon_, gcd(), gcdv(), CoinPackedVector::getElements(), CoinPackedVector::getIndices(), CoinPackedVector::getNumElements(), CoinPackedVector::insert(), power10ToMakeDoubleAnInt(), and CoinPackedVector::setVector().
00022 { 00023 int nRows=si.getNumRows(); // number of rows in the coefficient matrix 00024 int nCols=si.getNumCols(); // number of columns in the coefficient matrix 00025 int rowIndex; // index into the constraint matrix stored in row 00026 // order 00027 CoinPackedVector irow; // "integer row": working space to hold the integer 00028 // <= inequality derived from the rowIndex-th 00029 // constraint 00030 double b=0; // working space for the rhs of integer <= inequality 00031 bool * negative= new bool[nCols]; // negative[i]= true if coefficient of the 00032 // ith variable is negative and false 00033 // otherwise 00034 int k; // dummy iterator variable 00035 for ( k=0; k<nCols; k++ ) negative[k] = false; 00036 00037 const CoinPackedMatrix * rowCopy = 00038 si.getMatrixByRow(); // row copy: matrix stored in row order 00039 00041 // Main loop: // 00042 // For every row in the matrix, // 00043 // if we can derive a valid <= inequality in integer variables, then // 00044 // try to construct a simple rounding cut from the integer inequality. // 00045 // Add the resulting cut to the set of cuts. // 00047 00048 for (rowIndex=0; rowIndex<nRows; rowIndex++){ 00049 00050 // Only look at tight rows 00051 // double * pi=ekk_rowduals(model); 00052 // if (fabs(pi[row]) < epsilon_){ 00053 // continue; 00054 // } 00055 00056 // Try to derive an <= inequality in integer variables from the row 00057 // by netting out the continuous variables. 00058 // Store the value and the sign of the coefficients separately: 00059 // irow.getElements() contains the absolute values of the coefficients. 00060 // negative is a boolean vector indicating the sign of the coeffcients. 00061 // b is the rhs of the <= integer inequality 00062 00063 if (!deriveAnIntegerRow( si, 00064 rowIndex, 00065 rowCopy->getVector(rowIndex), 00066 irow, b, negative)) 00067 { 00068 00069 // Reset local data for the next iteration of the rowIndex-loop 00070 for(k=0; k<irow.getNumElements(); k++) negative[irow.getIndices()[k]]=false; 00071 irow.setVector(0,NULL,NULL); 00072 continue; 00073 } 00074 00075 // Euclid's greatest common divisor (gcd) algorithm applies to positive 00076 // INTEGERS. 00077 // Determine the power of 10 needed, so that multipylying the integer 00078 // inequality through by 10**power makes all coefficients essentially 00079 // integral. 00080 int power = power10ToMakeDoubleAnInt(irow.getNumElements(),irow.getElements(),epsilon_*1.0e-4); 00081 00082 // Now a vector to store the integer-ized values. For instance, 00083 // if x[i] is .66 and power is 1000 then xInt[i] will be 660 00084 int * xInt = NULL; 00085 if (power >=0) { 00086 00087 xInt = new int[irow.getNumElements()]; 00088 double dxInt; // a double version of xInt for error trapping 00089 00090 00091 #ifdef CGL_DEBUG 00092 printf("The (double) coefficients and their integer-ized counterparts:\n"); 00093 #endif 00094 00095 for (k=0; k<irow.getNumElements(); k++){ 00096 dxInt = irow.getElements()[k]*pow(10.0,power); 00097 xInt[k]= (int) (dxInt+0.5); // Need to add the 0.5 00098 // so that a dxInt=9.999 will give a xInt=1 00099 00100 #ifdef CGL_DEBUG 00101 printf("%g %g \n",irow.getElements()[k],dxInt); 00102 #endif 00103 00104 } 00105 00106 } else { 00107 00108 // If overflow is detected, one warning message is printed and 00109 // the row is skipped. 00110 #ifdef CGL_DEBUG 00111 printf("SimpleRounding: Warning: Overflow detected \n"); 00112 printf(" on %i of vars in processing row %i. Row skipped.\n", 00113 -power, rowIndex); 00114 #endif 00115 // reset local data for next iteration 00116 for(k=0; k<irow.getNumElements(); k++) negative[irow.getIndices()[k]]=false; 00117 irow.setVector(0,NULL,NULL); 00118 continue; 00119 } 00120 00121 // find greatest common divisor of the irow.elements 00122 int gcd = gcdv(irow.getNumElements(), xInt); 00123 00124 #ifdef CGL_DEBUG 00125 printf("The gcd of xInt is %i\n",gcd); 00126 #endif 00127 00128 // construct new cut by dividing through by gcd and 00129 // rounding down rhs and accounting for negatives 00130 CoinPackedVector cut; 00131 for (k=0; k<irow.getNumElements(); k++){ 00132 cut.insert(irow.getIndices()[k],xInt[k]/gcd); 00133 } 00134 double cutRhs = floor((b*pow(10.0,power))/gcd); 00135 00136 // un-negate the negated variables in the cut 00137 { 00138 const int s = cut.getNumElements(); 00139 const int * indices = cut.getIndices(); 00140 double* elements = cut.getElements(); 00141 for (k=0; k<s; k++){ 00142 int column=indices[k]; 00143 if (negative[column]) { 00144 elements[k] *= -1; 00145 } 00146 } 00147 } 00148 00149 // Create the row cut and add it to the set of cuts 00150 // It may not be violated 00151 if (fabs(cutRhs*gcd-b)> epsilon_){ // if the cut and row are different. 00152 OsiRowCut rc; 00153 rc.setRow(cut.getNumElements(),cut.getIndices(),cut.getElements()); 00154 rc.setLb(-DBL_MAX); 00155 rc.setUb(cutRhs); 00156 cs.insert(rc); 00157 00158 #ifdef CGL_DEBUG 00159 printf("Row %i had a simple rounding cut:\n",rowIndex); 00160 printf("Cut size: %i Cut rhs: %g Index Element \n", 00161 cut.getNumElements(), cutRhs); 00162 for (k=0; k<cut.getNumElements(); k++){ 00163 printf("%i %g\n",cut.getIndices()[k], cut.getElements()[k]); 00164 } 00165 printf("\n"); 00166 #endif 00167 } 00168 00169 // Reset local data for the next iteration of the rowIndex-loop 00170 for(k=0; k<irow.getNumElements(); k++) negative[irow.getIndices()[k]]=false; 00171 irow.setVector(0,NULL,NULL); 00172 delete [] xInt; 00173 00174 00175 } 00176 00177 delete [] negative; 00178 } |
|
Given a vector of doubles, x, with size elements and a positive tolerance, dataTol, this method returns the smallest power of 10 needed so that x[i]*10**power "is integer" for all i=0,...,size-1. change of definition of dataTol so that it refers to original data, not to scaled data as that seems to lead to problems. So if xScaled is x[i]*10**power and xInt is rounded(xScaled) then fabs(xScaled-xInt) <= dataTol*10**power. This means that dataTol should be smaller - say 1.0e-12 rather tahn 1.0e-8 Returns -number of times overflowed if the power is so big that it will cause overflow (i.e. integer stored will be bigger than 2**31). Test in cut generator. Definition at line 330 of file CglSimpleRounding.cpp. Referenced by generateCuts().
00335 { 00336 // Assumption: data precision is positive 00337 assert( dataTol > 0 ); 00338 00339 00340 int i; // loop iterator 00341 int maxPower=0; // maximum power of 10 used to convert any x[i] to an 00342 // integer 00343 // this is the number we are after. 00344 int power = 0; // power of 10 used to convert a particular x[i] to an 00345 // integer 00346 00347 #ifdef OLD_MULT 00348 double intPart; // the integer part of the number 00349 #endif 00350 double fracPart; // the fractional part of the number 00351 // we keep multiplying by 10 until the fractional part is 0 00352 // (well, really just until the factional part is less than 00353 // dataTol) 00354 00355 // JJF - code seems to fail sometimes as multiplying by 10 - so 00356 // definition of dataTol changed - see header file 00357 00358 const double multiplier[16]={1.0,1.0e1,1.0e2,1.0e3,1.0e4,1.0e5, 00359 1.0e6,1.0e7,1.0e8,1.0e9,1.0e10,1.0e11, 00360 1.0e12,1.0e13,1.0e14,1.0e15}; 00361 00362 // Loop through every element in the array in x 00363 for (i=0; i<size; i++){ 00364 power = 0; 00365 00366 #ifdef OLD_MULT 00367 // look at the fractional part of x[i] 00368 // FYI: if you want to modify this member function to take an input 00369 // vector x of arbitary sign, change this line below to 00370 // fracPart = modf(fabs(x[i]),&intPart); 00371 fracPart = modf(x[i],&intPart); 00372 00373 // if the fractional part is close enough to 0 or 1, we're done with this 00374 // value 00375 while(!(fracPart < dataTol || 1-fracPart < dataTol )) { 00376 // otherwise, multiply by 10 and look at the fractional part of the 00377 // result. 00378 ++power; 00379 fracPart = fracPart*10.0; 00380 fracPart = modf(fracPart,&intPart); 00381 } 00382 #else 00383 // use fabs as safer and does no harm 00384 double value = fabs(x[i]); 00385 double scaledValue; 00386 // Do loop - always using original value to stop round off error. 00387 // If we don't find in 15 goes give up 00388 for (power=0;power<16;power++) { 00389 double tolerance = dataTol*multiplier[power]; 00390 scaledValue = value*multiplier[power]; 00391 fracPart = scaledValue-floor(scaledValue); 00392 if(fracPart < tolerance || 1.0-fracPart < tolerance ) { 00393 break; 00394 } 00395 } 00396 if (power==16||scaledValue>2147483647) { 00397 #ifdef CGL_DEBUG 00398 printf("Overflow %g => %g, power %d\n",x[i],scaledValue,power); 00399 #endif 00400 return -1; 00401 } 00402 #endif 00403 #ifdef CGL_DEBUG 00404 printf("The smallest power of 10 to make %g integral = %i\n",x[i],power); 00405 #endif 00406 00407 00408 // keep track of the largest power needed so that at the end of the for 00409 // loop 00410 // x[i]*10**maxPower will be integral for all i 00411 if (maxPower < power) maxPower=power; 00412 } 00413 00414 return maxPower; 00415 } |
|
A function that tests the methods in the CglSimpleRounding class. The only reason for it not to be a member method is that this way it doesn't have to be compiled into the library. And that's a gain, because the library should be compiled with optimization on, but this method should be compiled with debugging. Definition at line 20 of file CglSimpleRoundingTest.cpp.
00023 { 00024 00025 // Test default constructor 00026 { 00027 CglSimpleRounding cg; 00028 } 00029 00030 // Test copy & assignment 00031 { 00032 CglSimpleRounding rhs; 00033 { 00034 CglSimpleRounding cg; 00035 CglSimpleRounding cgC(cg); 00036 rhs=cg; 00037 } 00038 } 00039 00040 // Test gcd and gcdn 00041 { 00042 CglSimpleRounding cg; 00043 int v = cg.gcd(122,356); 00044 assert(v==2); 00045 v=cg.gcd(356,122); 00046 assert(v==2); 00047 v=cg.gcd(54,67); 00048 assert(v==1); 00049 v=cg.gcd(67,54); 00050 assert(v==1); 00051 v=cg.gcd(485,485); 00052 assert(v==485); 00053 v=cg.gcd(17*13,17*23); 00054 assert( v==17); 00055 v=cg.gcd(17*13*5,17*23); 00056 assert( v==17); 00057 v=cg.gcd(17*13*23,17*23); 00058 assert(v==17*23); 00059 00060 int a[4] = {12, 20, 32, 400}; 00061 v= cg.gcdv(4,a); 00062 assert(v== 4); 00063 int b[4] = {782, 4692, 51, 2754}; 00064 v= cg.gcdv(4,b); 00065 assert(v== 17); 00066 int c[4] = {50, 40, 30, 10}; 00067 v= cg.gcdv(4,c); 00068 assert(v== 10); 00069 } 00070 00071 00072 // Test generate cuts method on exmip1.5.mps 00073 { 00074 CglSimpleRounding cg; 00075 00076 OsiSolverInterface * siP = baseSiP->clone(); 00077 std::string fn = mpsDir+"exmip1.5"; 00078 siP->readMps(fn.c_str(),"mps"); 00079 OsiCuts cuts; 00080 cg.generateCuts(*siP,cuts); 00081 00082 // there should be 3 cuts 00083 int nRowCuts = cuts.sizeRowCuts(); 00084 assert(nRowCuts==3); 00085 00086 // get the last "sr"=simple rounding cut that was derived 00087 OsiRowCut srRowCut2 = cuts.rowCut(2); 00088 CoinPackedVector srRowCutPV2 = srRowCut2.row(); 00089 00090 // this is what the last cut should look like: i.e. the "solution" 00091 const int solSize = 2; 00092 int solCols[solSize]={2,3}; 00093 double solCoefs[solSize]={5.0, 4.0}; 00094 OsiRowCut solRowCut; 00095 solRowCut.setRow(solSize,solCols,solCoefs); 00096 solRowCut.setLb(-DBL_MAX); 00097 solRowCut.setUb(2.0); 00098 00099 // Test for equality between the derived cut and the solution cut 00100 00101 // Note: testing two OsiRowCuts are equal invokes testing two 00102 // CoinPackedVectors are equal which invokes testing two doubles 00103 // are equal. Usually not a good idea to test that two doubles are equal, 00104 // but in this cut the "doubles" represent integer values. 00105 assert(srRowCut2 == solRowCut); 00106 00107 delete siP; 00108 } 00109 00110 // Test generate cuts method on p0033 00111 { 00112 CglSimpleRounding cg; 00113 00114 OsiSolverInterface * siP = baseSiP->clone(); 00115 std::string fn = mpsDir+"p0033"; 00116 siP->readMps(fn.c_str(),"mps"); 00117 OsiCuts cuts; 00118 cg.generateCuts(*siP,cuts); 00119 00120 // p0033 is the optimal solution to p0033 00121 int objIndices[14] = { 00122 0, 6, 7, 9, 13, 17, 18, 00123 22, 24, 25, 26, 27, 28, 29 }; 00124 CoinPackedVector p0033(14,objIndices,1.0); 00125 00126 // test that none of the generated cuts 00127 // chops off the optimal solution 00128 int nRowCuts = cuts.sizeRowCuts(); 00129 OsiRowCut rcut; 00130 CoinPackedVector rpv; 00131 int i; 00132 for (i=0; i<nRowCuts; i++){ 00133 rcut = cuts.rowCut(i); 00134 rpv = rcut.row(); 00135 double p0033Sum = (rpv*p0033).sum(); 00136 double rcutub = rcut.ub(); 00137 assert (p0033Sum <= rcutub); 00138 } 00139 00140 // test that the cuts improve the 00141 // lp objective function value 00142 siP->initialSolve(); 00143 double lpRelaxBefore=siP->getObjValue(); 00144 OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts); 00145 siP->resolve(); 00146 double lpRelaxAfter=siP->getObjValue(); 00147 #ifdef CGL_DEBUG 00148 printf("\n\nOrig LP min=%f\n",lpRelaxBefore); 00149 printf("Final LP min=%f\n\n",lpRelaxAfter); 00150 #endif 00151 assert( lpRelaxBefore < lpRelaxAfter ); 00152 00153 delete siP; 00154 00155 } 00156 00157 00158 } |