[1770] | 1 | #include <cstdlib> |
---|
| 2 | #include <cmath> |
---|
| 3 | #include <ctime> |
---|
| 4 | #include <cassert> |
---|
| 5 | #include <cstdio> |
---|
| 6 | #include <cstring> |
---|
| 7 | #include <algorithm> |
---|
| 8 | #include <vector> |
---|
| 9 | #include <string> |
---|
| 10 | #include <map> |
---|
| 11 | #include <OsiSolverInterface.hpp> |
---|
| 12 | #include "CbcMessage.hpp" |
---|
[1943] | 13 | #include "CbcHeuristic.hpp" |
---|
[1770] | 14 | #include <CbcModel.hpp> |
---|
| 15 | #include "CbcMipStartIO.hpp" |
---|
| 16 | #include "CoinTime.hpp" |
---|
| 17 | |
---|
| 18 | using namespace std; |
---|
| 19 | |
---|
| 20 | |
---|
| 21 | bool isNumericStr( const char *str ) |
---|
| 22 | { |
---|
[1886] | 23 | const size_t l = strlen(str); |
---|
[1770] | 24 | |
---|
[1886] | 25 | for ( size_t i=0 ; i<l ; ++i ) |
---|
[1943] | 26 | if (!(isdigit(str[i])||(str[i]=='.')||(str[i]=='-')||(str[i]=='e'))) |
---|
[1770] | 27 | return false; |
---|
| 28 | |
---|
| 29 | return true; |
---|
| 30 | } |
---|
| 31 | |
---|
| 32 | int readMIPStart( CbcModel * model, const char *fileName, |
---|
| 33 | vector< pair< string, double > > &colValues, |
---|
| 34 | double &/*solObj*/ ) |
---|
| 35 | { |
---|
| 36 | #define STR_SIZE 256 |
---|
| 37 | FILE *f = fopen( fileName, "r" ); |
---|
| 38 | if (!f) |
---|
| 39 | return 1; |
---|
| 40 | char line[STR_SIZE]; |
---|
| 41 | |
---|
| 42 | int nLine = 0; |
---|
| 43 | char printLine[STR_SIZE]; |
---|
| 44 | while (fgets( line, STR_SIZE, f )) |
---|
| 45 | { |
---|
| 46 | ++nLine; |
---|
| 47 | char col[4][STR_SIZE]; |
---|
| 48 | int nread = sscanf( line, "%s %s %s %s", col[0], col[1], col[2], col[3] ); |
---|
| 49 | if (!nread) |
---|
| 50 | continue; |
---|
| 51 | /* line with variable value */ |
---|
| 52 | if (strlen(col[0])&&isdigit(col[0][0])&&(nread>=3)) |
---|
| 53 | { |
---|
| 54 | if (!isNumericStr(col[0])) |
---|
| 55 | { |
---|
| 56 | sprintf( printLine, "Reading: %s, line %d - first column in mipstart file should be numeric, ignoring.", fileName, nLine ); |
---|
[2186] | 57 | model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol; |
---|
[1770] | 58 | continue; |
---|
| 59 | } |
---|
| 60 | if (!isNumericStr(col[2])) |
---|
| 61 | { |
---|
| 62 | sprintf( printLine, "Reading: %s, line %d - Third column in mipstart file should be numeric, ignoring.", fileName, nLine ); |
---|
[2186] | 63 | model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol; |
---|
[1770] | 64 | continue; |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | char *name = col[1]; |
---|
| 68 | double value = atof( col[2] ); |
---|
| 69 | |
---|
| 70 | colValues.push_back( pair<string, double>(string(name),value) ); |
---|
| 71 | } |
---|
| 72 | } |
---|
| 73 | |
---|
| 74 | if (colValues.size()) { |
---|
[2186] | 75 | sprintf( printLine,"MIPStart values read for %d variables.", static_cast<int>(colValues.size()) ); |
---|
| 76 | model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol; |
---|
[2013] | 77 | if (colValues.size()<model->getNumCols()) { |
---|
[2186] | 78 | int numberColumns = model->getNumCols(); |
---|
| 79 | OsiSolverInterface *solver = model->solver(); |
---|
| 80 | vector< pair< string, double > > fullValues; |
---|
| 81 | /* for fast search of column names */ |
---|
| 82 | map< string, int > colIdx; |
---|
| 83 | for (int i=0;i<numberColumns;i++) { |
---|
| 84 | fullValues.push_back( pair<string, double>(solver->getColName(i),0.0) ); |
---|
| 85 | colIdx[solver->getColName(i)] = i; |
---|
| 86 | } |
---|
| 87 | for ( int i=0 ; (i<static_cast<int>(colValues.size())) ; ++i ) { |
---|
| 88 | map< string, int >::const_iterator mIt = colIdx.find( colValues[i].first ); |
---|
| 89 | if ( mIt != colIdx.end() ) { |
---|
| 90 | const int idx = mIt->second; |
---|
| 91 | double v = colValues[i].second; |
---|
| 92 | fullValues[idx].second=v; |
---|
| 93 | } |
---|
| 94 | } |
---|
| 95 | colValues=fullValues; |
---|
[2013] | 96 | } |
---|
[2186] | 97 | } |
---|
| 98 | else { |
---|
[1770] | 99 | sprintf( printLine, "No mipstart solution read from %s", fileName ); |
---|
[2186] | 100 | model->messageHandler()->message(CBC_GENERAL, model->messages()) << printLine << CoinMessageEol; |
---|
[1770] | 101 | return 1; |
---|
| 102 | } |
---|
| 103 | |
---|
| 104 | fclose(f); |
---|
| 105 | return 0; |
---|
| 106 | } |
---|
| 107 | |
---|
| 108 | int computeCompleteSolution( CbcModel * model, |
---|
| 109 | const vector< string > colNames, |
---|
| 110 | const std::vector< std::pair< std::string, double > > &colValues, |
---|
| 111 | double *sol, double &obj ) |
---|
| 112 | { |
---|
[2187] | 113 | if (!model->getNumCols()) |
---|
| 114 | return 0; |
---|
| 115 | |
---|
[1770] | 116 | int status = 0; |
---|
| 117 | double compObj = COIN_DBL_MAX; |
---|
| 118 | bool foundIntegerSol = false; |
---|
| 119 | OsiSolverInterface *lp = model->solver()->clone(); |
---|
| 120 | map< string, int > colIdx; |
---|
[2186] | 121 | assert( (static_cast<int>(colNames.size())) == lp->getNumCols() ); |
---|
[1770] | 122 | /* for fast search of column names */ |
---|
[2186] | 123 | for ( int i=0 ; (i<static_cast<int>(colNames.size())) ; ++i ) |
---|
[1770] | 124 | colIdx[colNames[i]] = i; |
---|
| 125 | |
---|
| 126 | char printLine[STR_SIZE]; |
---|
| 127 | int fixed = 0; |
---|
| 128 | int notFound = 0; |
---|
| 129 | char colNotFound[256] = ""; |
---|
[1943] | 130 | int nContinuousFixed = 0; |
---|
[2186] | 131 | |
---|
[1943] | 132 | #ifndef JUST_FIX_INTEGER |
---|
| 133 | #define JUST_FIX_INTEGER 0 |
---|
| 134 | #endif |
---|
[2186] | 135 | |
---|
[1943] | 136 | #if JUST_FIX_INTEGER > 1 |
---|
| 137 | // all not mentioned are at zero |
---|
| 138 | for ( int i=0 ; (i<lp->getNumCols()) ; ++i ) |
---|
[2186] | 139 | { |
---|
[1943] | 140 | if (lp->isInteger(i)) |
---|
| 141 | lp->setColBounds( i, 0.0, 0.0 ); |
---|
[2186] | 142 | } |
---|
[1943] | 143 | #endif |
---|
[2186] | 144 | for ( int i=0 ; (i<static_cast<int>(colValues.size())) ; ++i ) |
---|
[1770] | 145 | { |
---|
| 146 | map< string, int >::const_iterator mIt = colIdx.find( colValues[i].first ); |
---|
| 147 | if ( mIt == colIdx.end() ) |
---|
| 148 | { |
---|
| 149 | if (!notFound) |
---|
| 150 | strcpy( colNotFound, colValues[i].first.c_str() ); |
---|
| 151 | notFound++; |
---|
| 152 | } |
---|
| 153 | else |
---|
| 154 | { |
---|
| 155 | const int idx = mIt->second; |
---|
| 156 | double v = colValues[i].second; |
---|
[1943] | 157 | #if JUST_FIX_INTEGER |
---|
| 158 | if (!lp->isInteger(idx)) |
---|
[2186] | 159 | continue; |
---|
[1943] | 160 | #endif |
---|
[1770] | 161 | if (v<1e-8) |
---|
| 162 | v = 0.0; |
---|
| 163 | if (lp->isInteger(idx)) // just to avoid small |
---|
| 164 | v = floor( v+0.5 ); // fractional garbage |
---|
[2186] | 165 | else |
---|
| 166 | nContinuousFixed++; |
---|
| 167 | |
---|
[1770] | 168 | lp->setColBounds( idx, v, v ); |
---|
| 169 | ++fixed; |
---|
| 170 | } |
---|
| 171 | } |
---|
| 172 | |
---|
| 173 | if (!fixed) |
---|
| 174 | { |
---|
| 175 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
[2186] | 176 | << "Warning: MIPstart solution is not valid, column names do not match, ignoring it." |
---|
| 177 | << CoinMessageEol; |
---|
[1770] | 178 | goto TERMINATE; |
---|
| 179 | } |
---|
| 180 | |
---|
[2186] | 181 | if ( notFound >= ( (static_cast<double>(colNames.size())) * 0.5 ) ) { |
---|
[1770] | 182 | sprintf( printLine, "Warning: %d column names were not found (e.g. %s) while filling solution.", notFound, colNotFound ); |
---|
[2186] | 183 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 184 | << printLine << CoinMessageEol; |
---|
[1770] | 185 | } |
---|
[1943] | 186 | #if JUST_FIX_INTEGER |
---|
| 187 | lp->setHintParam(OsiDoPresolveInInitial, true, OsiHintDo) ; |
---|
| 188 | #endif |
---|
[2013] | 189 | lp->setDblParam(OsiDualObjectiveLimit,COIN_DBL_MAX); |
---|
[1770] | 190 | lp->initialSolve(); |
---|
[2186] | 191 | |
---|
| 192 | if ( (lp->isProvenPrimalInfeasible()) || (lp->isProvenDualInfeasible()) ) |
---|
[1770] | 193 | { |
---|
[1943] | 194 | if (nContinuousFixed) { |
---|
[2186] | 195 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 196 | << "Trying just fixing integer variables." << CoinMessageEol; |
---|
| 197 | int numberColumns = lp->getNumCols(); |
---|
| 198 | const double *oldLower = model->solver()->getColLower(); |
---|
| 199 | const double *oldUpper = model->solver()->getColUpper(); |
---|
| 200 | for ( int i=0 ; i<numberColumns ; ++i ) { |
---|
| 201 | if (!lp->isInteger(i)) { |
---|
| 202 | lp->setColLower(i,oldLower[i]); |
---|
| 203 | lp->setColUpper(i,oldUpper[i]); |
---|
| 204 | } |
---|
| 205 | } |
---|
| 206 | |
---|
| 207 | lp->initialSolve(); |
---|
[1943] | 208 | } |
---|
[2186] | 209 | else |
---|
| 210 | { |
---|
| 211 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 212 | << "Fixing only non-zero variables." << CoinMessageEol; |
---|
| 213 | /* unfix all variables which are zero */ |
---|
| 214 | int notZeroAnymore = 0; |
---|
| 215 | for ( int i=0 ; (i<lp->getNumCols()) ; ++i ) |
---|
| 216 | if ( ((fabs(lp->getColLower()[i])) <= 1e-8) && (fabs(lp->getColLower()[i]-lp->getColUpper()[i]) <= 1e-8) ) |
---|
| 217 | { |
---|
| 218 | const double *oldLower = model->solver()->getColLower(); |
---|
| 219 | const double *oldUpper = model->solver()->getColUpper(); |
---|
| 220 | lp->setColLower(i,oldLower[i]); |
---|
| 221 | lp->setColUpper(i,oldUpper[i]); |
---|
| 222 | notZeroAnymore++; |
---|
| 223 | } |
---|
| 224 | if (notZeroAnymore) |
---|
| 225 | lp->initialSolve(); |
---|
[1943] | 226 | } |
---|
[1770] | 227 | } |
---|
| 228 | |
---|
[2186] | 229 | if (!lp->isProvenOptimal()) |
---|
| 230 | { |
---|
| 231 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 232 | << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol; |
---|
| 233 | status = 1; |
---|
| 234 | goto TERMINATE; |
---|
| 235 | } |
---|
| 236 | |
---|
[1770] | 237 | /* some additional effort is needed to provide an integer solution */ |
---|
| 238 | if ( lp->getFractionalIndices().size() > 0 ) |
---|
| 239 | { |
---|
[2186] | 240 | sprintf( printLine,"MIPStart solution provided values for %d of %d integer variables, %d variables are still fractional.", fixed, lp->getNumIntegers(), static_cast<int>(lp->getFractionalIndices().size()) ); |
---|
[1770] | 241 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 242 | << printLine << CoinMessageEol; |
---|
| 243 | double start = CoinCpuTime(); |
---|
[1943] | 244 | #if 1 |
---|
| 245 | CbcSerendipity heuristic(*model); |
---|
| 246 | heuristic.setFractionSmall(2.0); |
---|
| 247 | heuristic.setFeasibilityPumpOptions(1008013); |
---|
| 248 | int returnCode = heuristic.smallBranchAndBound(lp, |
---|
| 249 | 1000, sol, |
---|
| 250 | compObj, |
---|
| 251 | model->getCutoff(), |
---|
| 252 | "ReduceInMIPStart"); |
---|
| 253 | if ((returnCode&1) != 0) { |
---|
| 254 | sprintf( printLine,"Mini branch and bound defined values for remaining variables in %.2f seconds.", |
---|
| 255 | CoinCpuTime()-start); |
---|
| 256 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 257 | << printLine << CoinMessageEol; |
---|
| 258 | foundIntegerSol = true; |
---|
| 259 | obj = compObj; |
---|
| 260 | } |
---|
| 261 | #else |
---|
[1770] | 262 | CbcModel babModel( *lp ); |
---|
[2186] | 263 | lp->writeLp("lessFix"); |
---|
| 264 | babModel.setLogLevel( 2 ); |
---|
| 265 | babModel.setMaximumNodes( 1000 ); |
---|
[1770] | 266 | babModel.setMaximumSeconds( 60 ); |
---|
| 267 | babModel.branchAndBound(); |
---|
| 268 | if (babModel.bestSolution()) |
---|
| 269 | { |
---|
| 270 | sprintf( printLine,"Mini branch and bound defined values for remaining variables in %.2f seconds.", |
---|
| 271 | CoinCpuTime()-start); |
---|
| 272 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 273 | << printLine << CoinMessageEol; |
---|
| 274 | copy( babModel.bestSolution(), babModel.bestSolution()+babModel.getNumCols(), sol ); |
---|
| 275 | foundIntegerSol = true; |
---|
| 276 | obj = compObj = babModel.getObjValue(); |
---|
| 277 | } |
---|
[1943] | 278 | #endif |
---|
[1770] | 279 | else |
---|
| 280 | { |
---|
[2186] | 281 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
| 282 | << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol; |
---|
[1770] | 283 | status = 1; |
---|
| 284 | goto TERMINATE; |
---|
| 285 | } |
---|
| 286 | } |
---|
| 287 | else |
---|
| 288 | { |
---|
| 289 | foundIntegerSol = true; |
---|
| 290 | obj = compObj = lp->getObjValue(); |
---|
| 291 | copy( lp->getColSolution(), lp->getColSolution()+lp->getNumCols(), sol ); |
---|
| 292 | } |
---|
| 293 | |
---|
| 294 | if ( foundIntegerSol ) |
---|
| 295 | { |
---|
[2186] | 296 | sprintf( printLine,"MIPStart provided solution with cost %g", compObj); |
---|
[1770] | 297 | model->messageHandler()->message(CBC_GENERAL, model->messages()) |
---|
[2186] | 298 | << printLine << CoinMessageEol; |
---|
[1957] | 299 | #if 0 |
---|
[1943] | 300 | { |
---|
| 301 | int numberColumns=lp->getNumCols(); |
---|
| 302 | double largestInfeasibility = 0.0; |
---|
| 303 | double primalTolerance ; |
---|
| 304 | double offset; |
---|
| 305 | lp->getDblParam(OsiObjOffset, offset); |
---|
| 306 | lp->getDblParam(OsiPrimalTolerance, primalTolerance) ; |
---|
| 307 | const double *objective = lp->getObjCoefficients() ; |
---|
| 308 | const double * rowLower = lp->getRowLower() ; |
---|
| 309 | const double * rowUpper = lp->getRowUpper() ; |
---|
| 310 | const double * columnLower = lp->getColLower() ; |
---|
| 311 | const double * columnUpper = lp->getColUpper() ; |
---|
| 312 | int numberRows = lp->getNumRows() ; |
---|
| 313 | double *rowActivity = new double[numberRows] ; |
---|
| 314 | memset(rowActivity, 0, numberRows*sizeof(double)) ; |
---|
| 315 | double *rowSum = new double[numberRows] ; |
---|
| 316 | memset(rowSum, 0, numberRows*sizeof(double)) ; |
---|
| 317 | const double * element = lp->getMatrixByCol()->getElements(); |
---|
| 318 | const int * row = lp->getMatrixByCol()->getIndices(); |
---|
| 319 | const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts(); |
---|
| 320 | const int * columnLength = lp->getMatrixByCol()->getVectorLengths(); |
---|
| 321 | const CoinPackedMatrix * rowCopy = lp->getMatrixByRow(); |
---|
| 322 | const int * column = rowCopy->getIndices(); |
---|
| 323 | const int * rowLength = rowCopy->getVectorLengths(); |
---|
| 324 | const CoinBigIndex * rowStart = rowCopy->getVectorStarts(); |
---|
| 325 | const double * elementByRow = rowCopy->getElements(); |
---|
| 326 | double objValue=-offset; |
---|
| 327 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 328 | double value = sol[iColumn]; |
---|
| 329 | if (lp->isInteger(iColumn)) |
---|
| 330 | assert (fabs(value-floor(value+0.5))<1.0e-6); |
---|
| 331 | objValue += value*objective[iColumn]; |
---|
| 332 | if (value>columnUpper[iColumn]) { |
---|
| 333 | if (value-columnUpper[iColumn]>1.0e-8) |
---|
| 334 | printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]); |
---|
| 335 | value=columnUpper[iColumn]; |
---|
| 336 | } else if (value<columnLower[iColumn]) { |
---|
| 337 | if (value-columnLower[iColumn]<-1.0e-8) |
---|
| 338 | printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]); |
---|
| 339 | value=columnLower[iColumn]; |
---|
| 340 | } |
---|
| 341 | if (value) { |
---|
| 342 | CoinBigIndex start = columnStart[iColumn]; |
---|
| 343 | CoinBigIndex end = start + columnLength[iColumn]; |
---|
| 344 | for (CoinBigIndex j = start; j < end; j++) { |
---|
| 345 | int iRow = row[j]; |
---|
| 346 | if (fabs(value)<1.0e-6&&fabs(value*element[j])>1.0e-5) |
---|
| 347 | printf("Column %d row %d value %.8g element %g %s\n", |
---|
| 348 | iColumn,iRow,value,element[j],lp->isInteger(iColumn) ? "integer" : ""); |
---|
| 349 | rowActivity[iRow] += value * element[j]; |
---|
| 350 | rowSum[iRow] += fabs(value * element[j]); |
---|
| 351 | } |
---|
| 352 | } |
---|
| 353 | } |
---|
| 354 | for (int i = 0 ; i < numberRows ; i++) { |
---|
| 355 | #if 0 //def CLP_INVESTIGATE |
---|
| 356 | double inf; |
---|
| 357 | inf = rowLower[i] - rowActivity[i]; |
---|
| 358 | if (inf > primalTolerance) |
---|
| 359 | printf("Row %d inf %g sum %g %g <= %g <= %g\n", |
---|
| 360 | i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]); |
---|
| 361 | inf = rowActivity[i] - rowUpper[i]; |
---|
| 362 | if (inf > primalTolerance) |
---|
| 363 | printf("Row %d inf %g sum %g %g <= %g <= %g\n", |
---|
| 364 | i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]); |
---|
| 365 | #endif |
---|
| 366 | double infeasibility = CoinMax(rowActivity[i]-rowUpper[i], |
---|
| 367 | rowLower[i]-rowActivity[i]); |
---|
| 368 | // but allow for errors |
---|
| 369 | double factor = CoinMax(1.0,rowSum[i]*1.0e-3); |
---|
| 370 | if (infeasibility>largestInfeasibility*factor) { |
---|
| 371 | largestInfeasibility = infeasibility/factor; |
---|
| 372 | printf("Cinf of %g on row %d sum %g scaled %g\n", |
---|
| 373 | infeasibility,i,rowSum[i],largestInfeasibility); |
---|
| 374 | if (infeasibility>1.0e10) { |
---|
| 375 | for (CoinBigIndex j=rowStart[i]; |
---|
| 376 | j<rowStart[i]+rowLength[i];j++) { |
---|
| 377 | printf("col %d element %g\n", |
---|
| 378 | column[j],elementByRow[j]); |
---|
| 379 | } |
---|
| 380 | } |
---|
| 381 | } |
---|
| 382 | } |
---|
| 383 | delete [] rowActivity ; |
---|
| 384 | delete [] rowSum; |
---|
| 385 | if (largestInfeasibility > 10.0*primalTolerance) |
---|
| 386 | printf("Clargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue); |
---|
| 387 | else |
---|
| 388 | printf("Cfeasible (%g) - obj %g\n", largestInfeasibility,objValue); |
---|
| 389 | } |
---|
[1957] | 390 | #endif |
---|
[1770] | 391 | for ( int i=0 ; (i<lp->getNumCols()) ; ++i ) |
---|
| 392 | { |
---|
[1943] | 393 | #if 0 |
---|
[1770] | 394 | if (sol[i]<1e-8) |
---|
| 395 | sol[i] = 0.0; |
---|
| 396 | else |
---|
| 397 | if (lp->isInteger(i)) |
---|
| 398 | sol[i] = floor( sol[i]+0.5 ); |
---|
[1943] | 399 | #else |
---|
| 400 | if (lp->isInteger(i)) { |
---|
[1957] | 401 | //if (fabs(sol[i] - floor( sol[i]+0.5 ))>1.0e-8) |
---|
| 402 | //printf("bad sol for %d - %.12g\n",i,sol[i]); |
---|
[1943] | 403 | sol[i] = floor( sol[i]+0.5 ); |
---|
| 404 | } |
---|
| 405 | #endif |
---|
[1770] | 406 | } |
---|
[1957] | 407 | #if 0 |
---|
[1943] | 408 | { |
---|
| 409 | int numberColumns=lp->getNumCols(); |
---|
| 410 | double largestInfeasibility = 0.0; |
---|
| 411 | double primalTolerance ; |
---|
| 412 | double offset; |
---|
| 413 | lp->getDblParam(OsiObjOffset, offset); |
---|
| 414 | lp->getDblParam(OsiPrimalTolerance, primalTolerance) ; |
---|
| 415 | const double *objective = lp->getObjCoefficients() ; |
---|
| 416 | const double * rowLower = lp->getRowLower() ; |
---|
| 417 | const double * rowUpper = lp->getRowUpper() ; |
---|
| 418 | const double * columnLower = lp->getColLower() ; |
---|
| 419 | const double * columnUpper = lp->getColUpper() ; |
---|
| 420 | int numberRows = lp->getNumRows() ; |
---|
| 421 | double *rowActivity = new double[numberRows] ; |
---|
| 422 | memset(rowActivity, 0, numberRows*sizeof(double)) ; |
---|
| 423 | double *rowSum = new double[numberRows] ; |
---|
| 424 | memset(rowSum, 0, numberRows*sizeof(double)) ; |
---|
| 425 | const double * element = lp->getMatrixByCol()->getElements(); |
---|
| 426 | const int * row = lp->getMatrixByCol()->getIndices(); |
---|
| 427 | const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts(); |
---|
| 428 | const int * columnLength = lp->getMatrixByCol()->getVectorLengths(); |
---|
| 429 | const CoinPackedMatrix * rowCopy = lp->getMatrixByRow(); |
---|
| 430 | const int * column = rowCopy->getIndices(); |
---|
| 431 | const int * rowLength = rowCopy->getVectorLengths(); |
---|
| 432 | const CoinBigIndex * rowStart = rowCopy->getVectorStarts(); |
---|
| 433 | const double * elementByRow = rowCopy->getElements(); |
---|
| 434 | double objValue=-offset; |
---|
| 435 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 436 | double value = sol[iColumn]; |
---|
| 437 | if (lp->isInteger(iColumn)) |
---|
| 438 | assert (fabs(value-floor(value+0.5))<1.0e-6); |
---|
| 439 | objValue += value*objective[iColumn]; |
---|
| 440 | if (value>columnUpper[iColumn]) { |
---|
| 441 | if (value-columnUpper[iColumn]>1.0e-8) |
---|
| 442 | printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]); |
---|
| 443 | value=columnUpper[iColumn]; |
---|
| 444 | } else if (value<columnLower[iColumn]) { |
---|
| 445 | if (value-columnLower[iColumn]<-1.0e-8) |
---|
| 446 | printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]); |
---|
| 447 | value=columnLower[iColumn]; |
---|
| 448 | } |
---|
| 449 | if (value) { |
---|
| 450 | CoinBigIndex start = columnStart[iColumn]; |
---|
| 451 | CoinBigIndex end = start + columnLength[iColumn]; |
---|
| 452 | for (CoinBigIndex j = start; j < end; j++) { |
---|
| 453 | int iRow = row[j]; |
---|
| 454 | rowActivity[iRow] += value * element[j]; |
---|
| 455 | rowSum[iRow] += fabs(value * element[j]); |
---|
| 456 | } |
---|
| 457 | } |
---|
| 458 | } |
---|
| 459 | for (int i = 0 ; i < numberRows ; i++) { |
---|
| 460 | #if 0 //def CLP_INVESTIGATE |
---|
| 461 | double inf; |
---|
| 462 | inf = rowLower[i] - rowActivity[i]; |
---|
| 463 | if (inf > primalTolerance) |
---|
| 464 | printf("Row %d inf %g sum %g %g <= %g <= %g\n", |
---|
| 465 | i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]); |
---|
| 466 | inf = rowActivity[i] - rowUpper[i]; |
---|
| 467 | if (inf > primalTolerance) |
---|
| 468 | printf("Row %d inf %g sum %g %g <= %g <= %g\n", |
---|
| 469 | i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]); |
---|
| 470 | #endif |
---|
| 471 | double infeasibility = CoinMax(rowActivity[i]-rowUpper[i], |
---|
| 472 | rowLower[i]-rowActivity[i]); |
---|
| 473 | // but allow for errors |
---|
| 474 | double factor = CoinMax(1.0,rowSum[i]*1.0e-3); |
---|
| 475 | if (infeasibility>largestInfeasibility*factor) { |
---|
| 476 | largestInfeasibility = infeasibility/factor; |
---|
| 477 | printf("Dinf of %g on row %d sum %g scaled %g\n", |
---|
| 478 | infeasibility,i,rowSum[i],largestInfeasibility); |
---|
| 479 | if (infeasibility>1.0e10) { |
---|
| 480 | for (CoinBigIndex j=rowStart[i]; |
---|
| 481 | j<rowStart[i]+rowLength[i];j++) { |
---|
| 482 | printf("col %d element %g\n", |
---|
| 483 | column[j],elementByRow[j]); |
---|
| 484 | } |
---|
| 485 | } |
---|
| 486 | } |
---|
| 487 | } |
---|
| 488 | delete [] rowActivity ; |
---|
| 489 | delete [] rowSum; |
---|
| 490 | if (largestInfeasibility > 10.0*primalTolerance) |
---|
| 491 | printf("Dlargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue); |
---|
| 492 | else |
---|
| 493 | printf("Dfeasible (%g) - obj %g\n", largestInfeasibility,objValue); |
---|
| 494 | } |
---|
[1957] | 495 | #endif |
---|
[1943] | 496 | #if JUST_FIX_INTEGER |
---|
| 497 | const double * oldLower = model->solver()->getColLower(); |
---|
| 498 | const double * oldUpper = model->solver()->getColUpper(); |
---|
| 499 | const double * dj = lp->getReducedCost(); |
---|
| 500 | int nNaturalLB=0; |
---|
| 501 | int nMaybeLB=0; |
---|
| 502 | int nForcedLB=0; |
---|
| 503 | int nNaturalUB=0; |
---|
| 504 | int nMaybeUB=0; |
---|
| 505 | int nForcedUB=0; |
---|
| 506 | int nOther=0; |
---|
| 507 | for ( int i=0 ; i<lp->getNumCols() ; ++i ) { |
---|
| 508 | if (lp->isInteger(i)) { |
---|
| 509 | if (sol[i]==oldLower[i]) { |
---|
| 510 | if (dj[i]>1.0e-5) |
---|
| 511 | nNaturalLB++; |
---|
| 512 | else if (dj[i]<-1.0e-5) |
---|
| 513 | nForcedLB++; |
---|
| 514 | else |
---|
| 515 | nMaybeLB++; |
---|
| 516 | } else if (sol[i]==oldUpper[i]) { |
---|
| 517 | if (dj[i]<-1.0e-5) |
---|
| 518 | nNaturalUB++; |
---|
| 519 | else if (dj[i]>1.0e-5) |
---|
| 520 | nForcedUB++; |
---|
| 521 | else |
---|
| 522 | nMaybeUB++; |
---|
| 523 | } else { |
---|
| 524 | nOther++; |
---|
| 525 | } |
---|
| 526 | } |
---|
| 527 | } |
---|
| 528 | printf("%d other, LB %d natural, %d neutral, %d forced, UB %d natural, %d neutral, %d forced\n", |
---|
| 529 | nOther,nNaturalLB,nMaybeLB,nForcedLB, |
---|
| 530 | nNaturalUB,nMaybeUB,nForcedUB=0); |
---|
| 531 | #endif |
---|
[1770] | 532 | } |
---|
| 533 | |
---|
| 534 | TERMINATE: |
---|
| 535 | delete lp; |
---|
| 536 | return status; |
---|
| 537 | } |
---|
| 538 | #undef STR_SIZE |
---|