00001
00002
00003 #if defined(_MSC_VER)
00004
00005 # pragma warning(disable:4786)
00006 #endif
00007 #include <cstdlib>
00008 #include <cmath>
00009 #include <cstdio>
00010 #include <cfloat>
00011 #include <cassert>
00012
00013 #include "CglSimpleRounding.hpp"
00014 #include "CoinPackedVector.hpp"
00015 #include "CoinSort.hpp"
00016 #include "CoinPackedMatrix.hpp"
00017
00018
00019 void
00020 CglSimpleRounding::generateCuts(const OsiSolverInterface & si,
00021 OsiCuts & cs) const
00022 {
00023 int nRows=si.getNumRows();
00024 int nCols=si.getNumCols();
00025 int rowIndex;
00026
00027 CoinPackedVector irow;
00028
00029
00030 double b=0;
00031 bool * negative= new bool[nCols];
00032
00033
00034 int k;
00035 for ( k=0; k<nCols; k++ ) negative[k] = false;
00036
00037 const CoinPackedMatrix * rowCopy =
00038 si.getMatrixByRow();
00039
00041
00042
00043
00044
00045
00047
00048 for (rowIndex=0; rowIndex<nRows; rowIndex++){
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 if (!deriveAnIntegerRow( si,
00064 rowIndex,
00065 rowCopy->getVector(rowIndex),
00066 irow, b, negative))
00067 {
00068
00069
00070 for(k=0; k<irow.getNumElements(); k++) negative[irow.getIndices()[k]]=false;
00071 irow.setVector(0,NULL,NULL);
00072 continue;
00073 }
00074
00075
00076
00077
00078
00079
00080 int power = power10ToMakeDoubleAnInt(irow.getNumElements(),irow.getElements(),epsilon_*1.0e-4);
00081
00082
00083
00084 int * xInt = NULL;
00085 if (power >=0) {
00086
00087 xInt = new int[irow.getNumElements()];
00088 double dxInt;
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);
00098
00099
00100 #ifdef CGL_DEBUG
00101 printf("%g %g \n",irow.getElements()[k],dxInt);
00102 #endif
00103
00104 }
00105
00106 } else {
00107
00108
00109
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
00116 for(k=0; k<irow.getNumElements(); k++) negative[irow.getIndices()[k]]=false;
00117 irow.setVector(0,NULL,NULL);
00118 continue;
00119 }
00120
00121
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
00129
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
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
00150
00151 if (fabs(cutRhs*gcd-b)> epsilon_){
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
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 }
00179
00180
00181
00182
00183
00184
00185
00186
00187 bool
00188 CglSimpleRounding::deriveAnIntegerRow(
00189 const OsiSolverInterface & si,
00190 int rowIndex,
00191 const CoinShallowPackedVector & matrixRow,
00192 CoinPackedVector & irow,
00193 double & b,
00194 bool * negative) const
00195 {
00196 irow.clear();
00197 int i;
00198 double sign=1.0;
00199
00200
00201 int sizeOfRow=matrixRow.getNumElements();
00202
00203
00204 const char rowsense = si.getRowSense()[rowIndex];
00205
00206
00207 if (rowsense=='E' || rowsense=='N') {
00208 return 0;
00209 }
00210
00211 if (rowsense=='L'){
00212 b=si.getRightHandSide()[rowIndex];
00213 }
00214
00215
00216 if (rowsense=='G'){
00217 b=-si.getRightHandSide()[rowIndex];
00218 sign=-1.0;
00219 }
00220
00221
00222
00223
00224
00225 if (rowsense=='R') {
00226 b=si.getRightHandSide()[rowIndex];
00227 }
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 const double * colupper = si.getColUpper();
00248 const double * collower = si.getColLower();
00249
00250 for (i=0; i<sizeOfRow; i++){
00251
00252 if ( !si.isInteger( matrixRow.getIndices()[i] ) ) {
00253
00254 if((sign*matrixRow.getElements()[i])<-epsilon_){
00255
00256 if (colupper[matrixRow.getIndices()[i]] < si.getInfinity()){
00257
00258 b=b-(sign*matrixRow.getElements()[i]*colupper[matrixRow.getIndices()[i]]);
00259 }
00260 else
00261 return 0;
00262 }
00263
00264 else if((sign*matrixRow.getElements()[i])>epsilon_){
00265
00266 if (collower[matrixRow.getIndices()[i]] > -si.getInfinity()){
00267
00268 b=b-(sign*matrixRow.getElements()[i]*collower[matrixRow.getIndices()[i]]);
00269 }
00270 else
00271 return 0;
00272 }
00273
00274
00275 }
00276
00277 else{
00278
00279 if (colupper[matrixRow.getIndices()[i]]- collower[matrixRow.getIndices()[i]]<
00280 epsilon_){
00281 b=b-(sign*matrixRow.getElements()[i]*colupper[matrixRow.getIndices()[i]]);
00282 }
00283
00284
00285 else {
00286 irow.insert(matrixRow.getIndices()[i],sign*matrixRow.getElements()[i]);
00287 }
00288 }
00289 }
00290
00291
00292 if(irow.getNumElements() == 0){
00293 return 0;
00294 }
00295
00296
00297
00298
00299
00300
00301 {
00302 const int s = irow.getNumElements();
00303 const int * indices = irow.getIndices();
00304 double * elements = irow.getElements();
00305 for(i=0; i<s; i++){
00306 if (elements[i] < -epsilon_) {
00307 negative[indices[i]]= true;
00308 elements[i] *= -1;
00309 }
00310 }
00311 }
00312
00313 return 1;
00314 }
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329 int
00330 CglSimpleRounding::power10ToMakeDoubleAnInt(
00331 int size,
00332 const double * x,
00333 double dataTol) const
00334
00335 {
00336
00337 assert( dataTol > 0 );
00338
00339
00340 int i;
00341 int maxPower=0;
00342
00343
00344 int power = 0;
00345
00346
00347 #ifdef OLD_MULT
00348 double intPart;
00349 #endif
00350 double fracPart;
00351
00352
00353
00354
00355
00356
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
00363 for (i=0; i<size; i++){
00364 power = 0;
00365
00366 #ifdef OLD_MULT
00367
00368
00369
00370
00371 fracPart = modf(x[i],&intPart);
00372
00373
00374
00375 while(!(fracPart < dataTol || 1-fracPart < dataTol )) {
00376
00377
00378 ++power;
00379 fracPart = fracPart*10.0;
00380 fracPart = modf(fracPart,&intPart);
00381 }
00382 #else
00383
00384 double value = fabs(x[i]);
00385 double scaledValue;
00386
00387
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
00409
00410
00411 if (maxPower < power) maxPower=power;
00412 }
00413
00414 return maxPower;
00415 }
00416
00417
00418
00419
00420 CglSimpleRounding::CglSimpleRounding ()
00421 :
00422 CglCutGenerator(),
00423 epsilon_(1.0e-08)
00424 {
00425
00426 }
00427
00428
00429
00430 CglSimpleRounding::CglSimpleRounding (
00431 const CglSimpleRounding & source)
00432 :
00433 CglCutGenerator(source),
00434 epsilon_(source.epsilon_)
00435 {
00436
00437 }
00438
00439
00440
00441
00442
00443 CglSimpleRounding::~CglSimpleRounding ()
00444 {
00445
00446 }
00447
00448
00449
00450
00451 CglSimpleRounding &
00452 CglSimpleRounding::operator=(
00453 const CglSimpleRounding& rhs)
00454 {
00455 if (this != &rhs) {
00456 CglCutGenerator::operator=(rhs);
00457 epsilon_=rhs.epsilon_;
00458 }
00459 return *this;
00460 }