[1271] | 1 | /* $Id: unitTestClp.cpp 1880 2013-04-01 17:09:22Z forrest $ */ |
---|
[19] | 2 | // Copyright (C) 2002, International Business Machines |
---|
| 3 | // Corporation and others. All Rights Reserved. |
---|
[1573] | 4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
[19] | 5 | |
---|
[1573] | 6 | |
---|
[19] | 7 | #include <cstdio> |
---|
| 8 | #include <string> |
---|
| 9 | #include <iostream> |
---|
[1622] | 10 | #include <iomanip> |
---|
[19] | 11 | |
---|
| 12 | #include "CoinTime.hpp" |
---|
[1668] | 13 | #include "CoinFileIO.hpp" |
---|
[274] | 14 | #include "CbcModel.hpp" |
---|
[1315] | 15 | #include "CbcHeuristic.hpp" |
---|
[641] | 16 | #include "CbcCutGenerator.hpp" |
---|
[1121] | 17 | #include "CbcBranchCut.hpp" |
---|
[1053] | 18 | #include "CglProbing.hpp" |
---|
[274] | 19 | #include "OsiClpSolverInterface.hpp" |
---|
[931] | 20 | #include "ClpFactorization.hpp" |
---|
[838] | 21 | #include "OsiRowCutDebugger.hpp" |
---|
[19] | 22 | //############################################################################# |
---|
| 23 | |
---|
| 24 | #ifdef NDEBUG |
---|
| 25 | #undef NDEBUG |
---|
| 26 | #endif |
---|
| 27 | |
---|
[765] | 28 | //############################################################################# |
---|
[19] | 29 | |
---|
| 30 | // Display message on stdout and stderr |
---|
| 31 | void testingMessage( const char * const msg ) |
---|
| 32 | { |
---|
[1286] | 33 | std::cout << msg; |
---|
| 34 | //cout <<endl <<"*****************************************" |
---|
| 35 | // <<endl <<msg <<endl; |
---|
[19] | 36 | } |
---|
| 37 | |
---|
[765] | 38 | //############################################################################# |
---|
[19] | 39 | |
---|
[768] | 40 | static inline bool CbcTestFile(const std::string name) |
---|
[274] | 41 | { |
---|
[1286] | 42 | FILE *fp = fopen(name.c_str(), "r"); |
---|
| 43 | if (fp) { |
---|
| 44 | fclose(fp); |
---|
| 45 | return true; |
---|
| 46 | } |
---|
| 47 | return false; |
---|
[768] | 48 | } |
---|
| 49 | |
---|
| 50 | //############################################################################# |
---|
| 51 | |
---|
| 52 | bool CbcTestMpsFile(std::string& fname) |
---|
| 53 | { |
---|
[1286] | 54 | if (CbcTestFile(fname)) { |
---|
| 55 | return true; |
---|
| 56 | } |
---|
| 57 | if (CbcTestFile(fname + ".mps")) { |
---|
| 58 | fname += ".mps"; |
---|
| 59 | return true; |
---|
| 60 | } |
---|
| 61 | if (CbcTestFile(fname + ".MPS")) { |
---|
| 62 | fname += ".MPS"; |
---|
| 63 | return true; |
---|
| 64 | } |
---|
[1668] | 65 | if (CoinFileInput::haveGzipSupport()) { |
---|
| 66 | if (CbcTestFile(fname + ".gz")) { |
---|
[1286] | 67 | return true; |
---|
[1668] | 68 | } |
---|
| 69 | if (CbcTestFile(fname + ".mps.gz")) { |
---|
[1286] | 70 | fname += ".mps"; |
---|
| 71 | return true; |
---|
[1668] | 72 | } |
---|
| 73 | if (CbcTestFile(fname + ".MPS.gz")) { |
---|
[1286] | 74 | fname += ".MPS"; |
---|
| 75 | return true; |
---|
[1668] | 76 | } |
---|
| 77 | if (CbcTestFile(fname + ".MPS.GZ")) { |
---|
| 78 | fname += ".MPS"; |
---|
| 79 | return true; |
---|
| 80 | } |
---|
[1286] | 81 | } |
---|
[1668] | 82 | if (CoinFileInput::haveBzip2Support()) { |
---|
| 83 | if (CbcTestFile(fname + ".bz2")) { |
---|
| 84 | return true; |
---|
| 85 | } |
---|
| 86 | if (CbcTestFile(fname + ".mps.bz2")) { |
---|
| 87 | fname += ".mps"; |
---|
| 88 | return true; |
---|
| 89 | } |
---|
| 90 | if (CbcTestFile(fname + ".MPS.bz2")) { |
---|
[1286] | 91 | fname += ".MPS"; |
---|
| 92 | return true; |
---|
[1668] | 93 | } |
---|
| 94 | if (CbcTestFile(fname + ".MPS.BZ2")) { |
---|
| 95 | fname += ".MPS"; |
---|
| 96 | return true; |
---|
| 97 | } |
---|
[1286] | 98 | } |
---|
| 99 | return false; |
---|
[768] | 100 | } |
---|
| 101 | |
---|
| 102 | //############################################################################# |
---|
[1622] | 103 | /* |
---|
| 104 | jjf: testSwitch -2 unitTest, -1 normal (==2) |
---|
| 105 | |
---|
| 106 | MiplibTest might be more appropriate. |
---|
| 107 | |
---|
| 108 | TestSwitch and stuff[6] together control how much of miplib is executed: |
---|
| 109 | For testSwitch set to: |
---|
| 110 | -2: solve p0033 and p0201 only (the unit test) |
---|
| 111 | -1: solve miplib sets #0 and #1 |
---|
| 112 | 0: solve nothing |
---|
| 113 | k: execute sets j:k, where j is determined by the value of stuff[6] |
---|
| 114 | The last parameter of PUSH_MPS specifies the test set membership. |
---|
| 115 | |
---|
| 116 | For -miplib, -extra2 sets testSwitch, -extra3 sets stuff[6]. The command |
---|
| 117 | cbc -extra2 -2 -miplib |
---|
| 118 | will execute the unit test on the miplib directory. |
---|
| 119 | |
---|
| 120 | dirMiplib should end in the directory separator character for the platform. |
---|
| 121 | |
---|
| 122 | If you want to activate the row cut debugger for a given problem, change the |
---|
| 123 | last parameter of the PUSH_MPS macro for the problem to true. |
---|
| 124 | |
---|
| 125 | Returns 0 if all goes well, -1 if the Miplib directory is missing, otherwise |
---|
| 126 | 100*(number with bad objective)+(number that exceeded node limit) |
---|
| 127 | */ |
---|
| 128 | int CbcClpUnitTest (const CbcModel &saveModel, const std::string &dirMiplib, |
---|
| 129 | int testSwitch, const double *stuff) |
---|
[768] | 130 | { |
---|
[1622] | 131 | // Stop Windows popup |
---|
| 132 | WindowsErrorPopupBlocker() ; |
---|
| 133 | unsigned int m ; |
---|
[768] | 134 | |
---|
[1622] | 135 | // Do an existence check. |
---|
| 136 | std::string test1 = dirMiplib+"p0033"; |
---|
| 137 | bool doTest = CbcTestMpsFile(test1); |
---|
| 138 | if (!doTest) { |
---|
| 139 | std::cout |
---|
| 140 | << "Not doing miplib run as can't find mps files." << std::endl |
---|
| 141 | << "Perhaps you're trying to read gzipped (.gz) files without libz?" |
---|
| 142 | << std::endl ; |
---|
| 143 | return (-1) ; |
---|
| 144 | } |
---|
[1645] | 145 | int dfltPrecision = static_cast<int>(std::cout.precision()) ; |
---|
[1622] | 146 | /* |
---|
| 147 | Set the range of problems to be tested. testSwitch = -2 is special and is |
---|
| 148 | picked up below. |
---|
| 149 | */ |
---|
| 150 | int loSet = 0 ; |
---|
| 151 | int hiSet = 0 ; |
---|
| 152 | if (testSwitch == -1) { |
---|
| 153 | loSet = 0 ; |
---|
| 154 | hiSet = 1 ; |
---|
| 155 | } else if (testSwitch >= 0) { |
---|
| 156 | loSet = static_cast<int>(stuff[6]) ; |
---|
| 157 | hiSet = testSwitch ; |
---|
| 158 | std::cout |
---|
| 159 | << "Solving miplib problems in sets " << loSet |
---|
| 160 | << ":" << hiSet << "." << std::endl ; |
---|
[1286] | 161 | } |
---|
[1622] | 162 | /* |
---|
| 163 | Vectors to hold test problem names and characteristics. |
---|
| 164 | */ |
---|
[1286] | 165 | std::vector<std::string> mpsName ; |
---|
| 166 | std::vector<int> nRows ; |
---|
| 167 | std::vector<int> nCols ; |
---|
| 168 | std::vector<double> objValueC ; |
---|
| 169 | std::vector<double> objValue ; |
---|
| 170 | std::vector<int> testSet ; |
---|
[1622] | 171 | std::vector<bool> rowCutDebugger ; |
---|
| 172 | /* |
---|
| 173 | A macro to make the vector creation marginally readable. Parameters are |
---|
| 174 | name, rows, columns, integer objective, continuous objective, set ID, |
---|
| 175 | row cut debugger |
---|
| 176 | |
---|
| 177 | To enable the row cut debugger for a given problem, change the last |
---|
| 178 | parameter to true. Don't forget to turn it off before committing changes! |
---|
| 179 | */ |
---|
[274] | 180 | #define PUSH_MPS(zz_mpsName_zz,\ |
---|
| 181 | zz_nRows_zz,zz_nCols_zz,zz_objValue_zz,zz_objValueC_zz, \ |
---|
[1622] | 182 | zz_testSet_zz, zz_rcDbg_zz) \ |
---|
[274] | 183 | mpsName.push_back(zz_mpsName_zz) ; \ |
---|
| 184 | nRows.push_back(zz_nRows_zz) ; \ |
---|
| 185 | nCols.push_back(zz_nCols_zz) ; \ |
---|
| 186 | objValueC.push_back(zz_objValueC_zz) ; \ |
---|
[1271] | 187 | testSet.push_back(zz_testSet_zz) ; \ |
---|
[1622] | 188 | objValue.push_back(zz_objValue_zz) ; \ |
---|
| 189 | rowCutDebugger.push_back(zz_rcDbg_zz) ; |
---|
| 190 | /* |
---|
| 191 | Push the miplib problems. Except for -2 (unitTest), push all, even if we're |
---|
| 192 | not going to do all of them. |
---|
| 193 | */ |
---|
| 194 | if (testSwitch == -2) { |
---|
| 195 | PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0, false); |
---|
| 196 | PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0, false); |
---|
| 197 | // PUSH_MPS("flugpl", 18, 18, 1201500, 1167185.7, 0, false); |
---|
| 198 | } else { |
---|
| 199 | /* |
---|
| 200 | Load up the problem vector. Note that the row counts here include the |
---|
| 201 | objective function. |
---|
| 202 | */ |
---|
| 203 | PUSH_MPS("10teams", 230, 2025, 924, 917, 1, false); |
---|
| 204 | PUSH_MPS("air03", 124, 10757, 340160, 338864.25, 0, false); |
---|
| 205 | PUSH_MPS("air04", 823, 8904, 56137, 55535.436, 2, false); |
---|
| 206 | PUSH_MPS("air05", 426, 7195, 26374, 25877.609, 2, false); |
---|
| 207 | PUSH_MPS("arki001", 1048, 1388, 7580813.0459, 7579599.80787, 7, false); |
---|
| 208 | PUSH_MPS("bell3a", 123, 133, 878430.32, 862578.64, 0, false); |
---|
| 209 | PUSH_MPS("bell5", 91, 104, 8966406.49, 8608417.95, 1, false); |
---|
| 210 | PUSH_MPS("blend2", 274, 353, 7.598985, 6.9156751140, 0, false); |
---|
| 211 | PUSH_MPS("cap6000", 2176, 6000, -2451377, -2451537.325, 1, false); |
---|
| 212 | PUSH_MPS("dano3mip", 3202, 13873, 728.1111, 576.23162474, 7, false); |
---|
| 213 | PUSH_MPS("danoint", 664, 521, 65.67, 62.637280418, 6, false); |
---|
| 214 | PUSH_MPS("dcmulti", 290, 548, 188182, 183975.5397, 0, false); |
---|
| 215 | PUSH_MPS("dsbmip", 1182, 1886, -305.19817501, -305.19817501, 0, false); |
---|
| 216 | PUSH_MPS("egout", 98, 141, 568.101, 149.589, 0, false); |
---|
| 217 | PUSH_MPS("enigma", 21, 100, 0.0, 0.0, 0, false); |
---|
| 218 | PUSH_MPS("fast0507", 507, 63009, 174, 172.14556668, 5, false); |
---|
| 219 | PUSH_MPS("fiber", 363, 1298, 405935.18000, 156082.51759, 0, false); |
---|
| 220 | PUSH_MPS("fixnet6", 478, 878, 3983, 1200.88, 1, false); |
---|
| 221 | PUSH_MPS("flugpl", 18, 18, 1201500, 1167185.7, 0, false); |
---|
| 222 | PUSH_MPS("gen", 780, 870, 112313, 112130.0, 0, false); |
---|
| 223 | PUSH_MPS("gesa2", 1392, 1224, 25779856.372, 25476489.678, 1, false); |
---|
| 224 | PUSH_MPS("gesa2_o", 1248, 1224, 25779856.372, 25476489.678, 1, false); |
---|
| 225 | PUSH_MPS("gesa3", 1368, 1152, 27991042.648, 27833632.451, 0, false); |
---|
| 226 | PUSH_MPS("gesa3_o", 1224, 1152, 27991042.648, 27833632.451, 0, false); |
---|
| 227 | PUSH_MPS("gt2", 29, 188, 21166.000, 13460.233074, 0, false); |
---|
| 228 | PUSH_MPS("harp2", 112, 2993, -73899798.00, -74353341.502, 6, false); |
---|
| 229 | PUSH_MPS("khb05250", 101, 1350, 106940226, 95919464.0, 0, false); |
---|
| 230 | PUSH_MPS("l152lav", 97, 1989, 4722, 4656.36, 1, false); |
---|
| 231 | PUSH_MPS("lseu", 28, 89, 1120, 834.68, 0, false); |
---|
| 232 | PUSH_MPS("mas74", 13, 151, 11801.18573, 10482.79528, 3, false); |
---|
| 233 | PUSH_MPS("mas76", 12, 151, 40005.05414, 38893.9036, 2, false); |
---|
| 234 | PUSH_MPS("misc03", 96, 160, 3360, 1910., 0, false); |
---|
| 235 | PUSH_MPS("misc06", 820, 1808, 12850.8607, 12841.6, 0, false); |
---|
| 236 | PUSH_MPS("misc07", 212, 260, 2810, 1415.0, 1, false); |
---|
| 237 | PUSH_MPS("mitre", 2054, 10724, 115155, 114740.5184, 1, false); |
---|
| 238 | PUSH_MPS("mkc", 3411, 5325, -553.75, -611.85, 7, false); // suboptimal |
---|
| 239 | PUSH_MPS("mod008", 6, 319, 307, 290.9, 0, false); |
---|
| 240 | PUSH_MPS("mod010", 146, 2655, 6548, 6532.08, 0, false); |
---|
| 241 | PUSH_MPS("mod011", 4480, 10958, -54558535, -62121982.55, 2, false); |
---|
| 242 | PUSH_MPS("modglob", 291, 422, 20740508, 20430947., 2, false); |
---|
| 243 | PUSH_MPS("noswot", 182, 128, -43, -43.0, 6, false); |
---|
| 244 | PUSH_MPS("nw04", 36, 87482, 16862, 16310.66667, 1, false); |
---|
| 245 | PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0, false); |
---|
| 246 | PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0, false); |
---|
| 247 | PUSH_MPS("p0282", 241, 282, 258411, 176867.50, 0, false); |
---|
| 248 | PUSH_MPS("p0548", 176, 548, 8691, 315.29, 0, false); |
---|
| 249 | PUSH_MPS("p2756", 755, 2756, 3124, 2688.75, 0, false); |
---|
| 250 | PUSH_MPS("pk1", 45, 86, 11.0, 0.0, 2, false); |
---|
| 251 | PUSH_MPS("pp08a", 136, 240, 7350.0, 2748.3452381, 1, false); |
---|
| 252 | PUSH_MPS("pp08aCUTS", 246, 240, 7350.0, 5480.6061563, 1, false); |
---|
| 253 | PUSH_MPS("qiu", 1192, 840, -132.873137, -931.638857, 3, false); |
---|
| 254 | PUSH_MPS("qnet1", 503, 1541, 16029.692681, 14274.102667, 0, false); |
---|
| 255 | PUSH_MPS("qnet1_o", 456, 1541, 16029.692681, 12095.571667, 0, false); |
---|
| 256 | PUSH_MPS("rentacar", 6803, 9557, 30356761, 28806137.644, 0, false); |
---|
| 257 | PUSH_MPS("rgn", 24, 180, 82.1999, 48.7999, 0, false); |
---|
| 258 | PUSH_MPS("rout", 291, 556, 1077.56, 981.86428571, 3, false); |
---|
| 259 | PUSH_MPS("set1ch", 492, 712, 54537.75, 32007.73, 5, false); |
---|
| 260 | PUSH_MPS("seymour", 4944, 1372, 423, 403.84647413, 7, false); |
---|
| 261 | PUSH_MPS("seymour_1", 4944, 1372, 410.76370, 403.84647413, 5, false); |
---|
| 262 | PUSH_MPS("stein27", 118, 27, 18, 13.0, 0, false); |
---|
| 263 | PUSH_MPS("stein45", 331, 45, 30, 22.0, 1, false); |
---|
| 264 | PUSH_MPS("swath", 884, 6805, 497.603, 334.4968581, 7, false); |
---|
| 265 | PUSH_MPS("vpm1", 234, 378, 20, 15.4167, 0, false); |
---|
| 266 | PUSH_MPS("vpm2", 234, 378, 13.75, 9.8892645972, 0, false); |
---|
| 267 | } |
---|
[274] | 268 | #undef PUSH_MPS |
---|
[1286] | 269 | |
---|
[1622] | 270 | /* |
---|
| 271 | Normally the problems are executed in order. Define RANDOM_ORDER below to |
---|
| 272 | randomize. |
---|
| 273 | |
---|
| 274 | #define RANDOM_ORDER |
---|
| 275 | */ |
---|
| 276 | int which[100]; |
---|
| 277 | int nLoop = static_cast<int>(mpsName.size()); |
---|
| 278 | assert (nLoop <= 100); |
---|
| 279 | for (int i = 0; i < nLoop; i++) which[i] = i; |
---|
| 280 | |
---|
| 281 | # ifdef RANDOM_ORDER |
---|
| 282 | unsigned int iTime = static_cast<unsigned int>(CoinGetTimeOfDay()-1.256e9); |
---|
| 283 | std::cout << "Time (seed) " << iTime << "." << std::endl ; |
---|
| 284 | double sort[100]; |
---|
| 285 | CoinDrand48(true,iTime); |
---|
| 286 | for (int i = 0; i < nLoop; i++) sort[i] = CoinDrand48(); |
---|
| 287 | CoinSort_2(sort,sort+nLoop,which); |
---|
| 288 | # endif |
---|
| 289 | |
---|
| 290 | int problemCnt = 0; |
---|
| 291 | for (m = 0 ; m < mpsName.size() ; m++) { |
---|
| 292 | int setID = testSet[m]; |
---|
| 293 | if (loSet <= setID && setID <= hiSet) problemCnt++; |
---|
| 294 | } |
---|
| 295 | |
---|
| 296 | int numberFailures = 0; |
---|
| 297 | int numberAttempts = 0; |
---|
| 298 | int numProbSolved = 0; |
---|
| 299 | double timeTaken = 0.0; |
---|
| 300 | |
---|
| 301 | //#define CLP_FACTORIZATION_INSTRUMENT |
---|
| 302 | # ifdef CLP_FACTORIZATION_INSTRUMENT |
---|
| 303 | double timeTakenFac = 0.0; |
---|
| 304 | # endif |
---|
| 305 | |
---|
| 306 | /* |
---|
| 307 | Open the main loop to step through the MPS problems. |
---|
| 308 | */ |
---|
| 309 | for (unsigned int mw = 0 ; mw < mpsName.size() ; mw++) { |
---|
| 310 | m = which[mw]; |
---|
| 311 | int setID = testSet[m]; |
---|
| 312 | // Skip if problem is not in specified problem set(s) |
---|
| 313 | if (!(loSet <= setID && setID <= hiSet)) continue ; |
---|
| 314 | |
---|
| 315 | numberAttempts++; |
---|
| 316 | std::cout << " processing mps file: " << mpsName[m] |
---|
| 317 | << " (" << numberAttempts << " out of " |
---|
| 318 | << problemCnt << ")" << std::endl ; |
---|
| 319 | /* |
---|
| 320 | Stage 1: Read the MPS and make sure the size of the constraint matrix |
---|
| 321 | is correct. |
---|
| 322 | */ |
---|
| 323 | CbcModel *model = new CbcModel(saveModel) ; |
---|
| 324 | |
---|
| 325 | std::string fn = dirMiplib+mpsName[m] ; |
---|
| 326 | if (!CbcTestMpsFile(fn)) { |
---|
| 327 | std::cout << "ERROR: Cannot find MPS file " << fn << "." << std::endl ; |
---|
| 328 | continue; |
---|
[1286] | 329 | } |
---|
[1622] | 330 | model->solver()->readMps(fn.c_str(), "") ; |
---|
| 331 | assert(model->getNumRows() == nRows[m]) ; |
---|
| 332 | assert(model->getNumCols() == nCols[m]) ; |
---|
[1286] | 333 | |
---|
[1622] | 334 | // Careful! We're initialising for the benefit of other code. |
---|
| 335 | CoinDrand48(true,1234567); |
---|
| 336 | //printf("RAND1 %g %g\n",CoinDrand48(true,1234567),model->randomNumberGenerator()->randomDouble()); |
---|
| 337 | //printf("RAND1 %g\n",CoinDrand48(true,1234567)); |
---|
[1286] | 338 | |
---|
[1622] | 339 | // Higher limits for the serious problems. |
---|
| 340 | int testMaximumNodes = 200000; |
---|
| 341 | if (hiSet > 1) |
---|
| 342 | testMaximumNodes = 20000000; |
---|
| 343 | if (model->getMaximumNodes() > testMaximumNodes) { |
---|
| 344 | model->setMaximumNodes(testMaximumNodes); |
---|
| 345 | } |
---|
| 346 | /* |
---|
| 347 | Stage 2: Call solver to solve the problem. |
---|
| 348 | */ |
---|
[1286] | 349 | |
---|
[1622] | 350 | # ifdef CLP_FACTORIZATION_INSTRUMENT |
---|
| 351 | extern double factorization_instrument(int type); |
---|
| 352 | double facTime1 = factorization_instrument(0); |
---|
| 353 | std::cout |
---|
| 354 | << "Factorization - initial solve " << facTime1 << " seconds." |
---|
| 355 | << std::endl ; |
---|
| 356 | timeTakenFac += facTime1; |
---|
| 357 | # endif |
---|
[1286] | 358 | |
---|
[1622] | 359 | double startTime = CoinCpuTime()+CoinCpuTimeJustChildren(); |
---|
[1286] | 360 | |
---|
[1622] | 361 | // Setup specific to clp |
---|
| 362 | OsiClpSolverInterface *siClp = |
---|
| 363 | dynamic_cast<OsiClpSolverInterface *>(model->solver()) ; |
---|
| 364 | ClpSimplex *modelC = NULL; |
---|
| 365 | if (siClp) { |
---|
| 366 | modelC = siClp->getModelPtr(); |
---|
| 367 | ClpMatrixBase * matrix = modelC->clpMatrix(); |
---|
| 368 | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
---|
| 369 | if (stuff && stuff[9] && clpMatrix) { |
---|
| 370 | // vector matrix! |
---|
| 371 | clpMatrix->makeSpecialColumnCopy(); |
---|
| 372 | } |
---|
[1286] | 373 | |
---|
[1622] | 374 | # ifdef JJF_ZERO |
---|
| 375 | if (clpMatrix) { |
---|
| 376 | int numberRows = clpMatrix->getNumRows(); |
---|
| 377 | int numberColumns = clpMatrix->getNumCols(); |
---|
| 378 | double * elements = clpMatrix->getMutableElements(); |
---|
| 379 | const int * row = clpMatrix->getIndices(); |
---|
| 380 | const CoinBigIndex * columnStart = clpMatrix->getVectorStarts(); |
---|
| 381 | const int * columnLength = clpMatrix->getVectorLengths(); |
---|
| 382 | double * smallest = new double [numberRows]; |
---|
| 383 | double * largest = new double [numberRows]; |
---|
| 384 | char * flag = new char [numberRows]; |
---|
| 385 | CoinZeroN(flag, numberRows); |
---|
| 386 | for (int i = 0; i < numberRows; i++) { |
---|
| 387 | smallest[i] = COIN_DBL_MAX; |
---|
| 388 | largest[i] = 0.0; |
---|
| 389 | } |
---|
| 390 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 391 | bool isInteger = modelC->isInteger(iColumn); |
---|
| 392 | CoinBigIndex j; |
---|
| 393 | for (j = columnStart[iColumn]; |
---|
| 394 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 395 | int iRow = row[j]; |
---|
| 396 | double value = fabs(elements[j]); |
---|
| 397 | if (!isInteger) |
---|
| 398 | flag[iRow] = 1; |
---|
| 399 | smallest[iRow] = CoinMin(smallest[iRow], value); |
---|
| 400 | largest[iRow] = CoinMax(largest[iRow], value); |
---|
| 401 | } |
---|
| 402 | } |
---|
| 403 | double * rowLower = modelC->rowLower(); |
---|
| 404 | double * rowUpper = modelC->rowUpper(); |
---|
| 405 | bool changed = false; |
---|
| 406 | for (int i = 0; i < numberRows; i++) { |
---|
| 407 | if (flag[i] && smallest[i] > 10.0 && false) { |
---|
| 408 | smallest[i] = 1.0 / smallest[i]; |
---|
| 409 | if (rowLower[i] > -1.0e20) |
---|
| 410 | rowLower[i] *= smallest[i]; |
---|
| 411 | if (rowUpper[i] < 1.0e20) |
---|
| 412 | rowUpper[i] *= smallest[i]; |
---|
| 413 | changed = true; |
---|
| 414 | } else { |
---|
| 415 | smallest[i] = 0.0; |
---|
| 416 | } |
---|
| 417 | } |
---|
| 418 | if (changed) { |
---|
| 419 | printf("SCALED\n"); |
---|
| 420 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 421 | CoinBigIndex j; |
---|
| 422 | for (j = columnStart[iColumn]; |
---|
| 423 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 424 | int iRow = row[j]; |
---|
| 425 | if (smallest[iRow]) |
---|
| 426 | elements[j] *= smallest[iRow]; |
---|
| 427 | } |
---|
| 428 | } |
---|
| 429 | } |
---|
| 430 | delete [] smallest; |
---|
| 431 | delete [] largest; |
---|
| 432 | delete [] flag; |
---|
| 433 | } |
---|
| 434 | # endif // JJF_ZERO |
---|
| 435 | |
---|
| 436 | model->checkModel(); |
---|
| 437 | modelC->tightenPrimalBounds(0.0, 0, true); |
---|
| 438 | model->initialSolve(); |
---|
| 439 | if (modelC->dualBound() == 1.0e10) { |
---|
| 440 | // user did not set - so modify |
---|
| 441 | // get largest scaled away from bound |
---|
| 442 | ClpSimplex temp = *modelC; |
---|
| 443 | temp.dual(0, 7); |
---|
| 444 | double largestScaled = 1.0e-12; |
---|
| 445 | double largest = 1.0e-12; |
---|
| 446 | int numberRows = temp.numberRows(); |
---|
| 447 | const double * rowPrimal = temp.primalRowSolution(); |
---|
| 448 | const double * rowLower = temp.rowLower(); |
---|
| 449 | const double * rowUpper = temp.rowUpper(); |
---|
| 450 | const double * rowScale = temp.rowScale(); |
---|
| 451 | int iRow; |
---|
| 452 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
| 453 | double value = rowPrimal[iRow]; |
---|
| 454 | double above = value - rowLower[iRow]; |
---|
| 455 | double below = rowUpper[iRow] - value; |
---|
| 456 | if (above < 1.0e12) { |
---|
| 457 | largest = CoinMax(largest, above); |
---|
| 458 | } |
---|
| 459 | if (below < 1.0e12) { |
---|
| 460 | largest = CoinMax(largest, below); |
---|
| 461 | } |
---|
| 462 | if (rowScale) { |
---|
| 463 | double multiplier = rowScale[iRow]; |
---|
| 464 | above *= multiplier; |
---|
| 465 | below *= multiplier; |
---|
| 466 | } |
---|
| 467 | if (above < 1.0e12) { |
---|
| 468 | largestScaled = CoinMax(largestScaled, above); |
---|
| 469 | } |
---|
| 470 | if (below < 1.0e12) { |
---|
| 471 | largestScaled = CoinMax(largestScaled, below); |
---|
| 472 | } |
---|
| 473 | } |
---|
| 474 | |
---|
| 475 | int numberColumns = temp.numberColumns(); |
---|
| 476 | const double * columnPrimal = temp.primalColumnSolution(); |
---|
| 477 | const double * columnLower = temp.columnLower(); |
---|
| 478 | const double * columnUpper = temp.columnUpper(); |
---|
| 479 | const double * columnScale = temp.columnScale(); |
---|
| 480 | int iColumn; |
---|
| 481 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 482 | double value = columnPrimal[iColumn]; |
---|
| 483 | double above = value - columnLower[iColumn]; |
---|
| 484 | double below = columnUpper[iColumn] - value; |
---|
| 485 | if (above < 1.0e12) { |
---|
| 486 | largest = CoinMax(largest, above); |
---|
| 487 | } |
---|
| 488 | if (below < 1.0e12) { |
---|
| 489 | largest = CoinMax(largest, below); |
---|
| 490 | } |
---|
| 491 | if (columnScale) { |
---|
| 492 | double multiplier = 1.0 / columnScale[iColumn]; |
---|
| 493 | above *= multiplier; |
---|
| 494 | below *= multiplier; |
---|
| 495 | } |
---|
| 496 | if (above < 1.0e12) { |
---|
| 497 | largestScaled = CoinMax(largestScaled, above); |
---|
| 498 | } |
---|
| 499 | if (below < 1.0e12) { |
---|
| 500 | largestScaled = CoinMax(largestScaled, below); |
---|
| 501 | } |
---|
| 502 | } |
---|
| 503 | std::cout << "Largest (scaled) away from bound " << largestScaled |
---|
| 504 | << " unscaled " << largest << std::endl; |
---|
| 505 | # ifdef JJF_ZERO |
---|
| 506 | modelC->setDualBound(CoinMax(1.0001e8, |
---|
| 507 | CoinMin(1000.0*largestScaled,1.00001e10))); |
---|
| 508 | # else |
---|
| 509 | modelC->setDualBound(CoinMax(1.0001e9, |
---|
| 510 | CoinMin(1000.0*largestScaled,1.0001e10))); |
---|
| 511 | # endif |
---|
| 512 | } |
---|
| 513 | } // end clp-specific setup |
---|
| 514 | /* |
---|
| 515 | Cut passes: For small models (n < 500) always do 100 passes, if possible |
---|
| 516 | (-100). For larger models, use minimum drop to stop (100, 20). |
---|
| 517 | */ |
---|
| 518 | model->setMinimumDrop(CoinMin(5.0e-2, |
---|
| 519 | fabs(model->getMinimizationObjValue())*1.0e-3+1.0e-4)); |
---|
| 520 | if (CoinAbs(model->getMaximumCutPassesAtRoot()) <= 100) { |
---|
| 521 | if (model->getNumCols() < 500) { |
---|
| 522 | model->setMaximumCutPassesAtRoot(-100); |
---|
| 523 | } else if (model->getNumCols() < 5000) { |
---|
| 524 | model->setMaximumCutPassesAtRoot(100); |
---|
| 525 | } else { |
---|
| 526 | model->setMaximumCutPassesAtRoot(20); |
---|
| 527 | } |
---|
| 528 | } |
---|
| 529 | // If defaults then increase trust for small models |
---|
| 530 | if (model->numberStrong() == 5 && model->numberBeforeTrust() == 10) { |
---|
| 531 | int numberColumns = model->getNumCols(); |
---|
| 532 | if (numberColumns <= 50) { |
---|
| 533 | model->setNumberBeforeTrust(1000); |
---|
| 534 | } else if (numberColumns <= 100) { |
---|
| 535 | model->setNumberBeforeTrust(100); |
---|
| 536 | } else if (numberColumns <= 300) { |
---|
| 537 | model->setNumberBeforeTrust(50); |
---|
| 538 | } |
---|
| 539 | } |
---|
| 540 | //if (model->getNumCols()>=500) { |
---|
| 541 | // switch off Clp stuff |
---|
| 542 | //model->setFastNodeDepth(-1); |
---|
| 543 | //} |
---|
| 544 | /* |
---|
| 545 | Activate the row cut debugger, if requested. |
---|
| 546 | */ |
---|
| 547 | if (rowCutDebugger[m] == true) { |
---|
| 548 | std::string probName ; |
---|
| 549 | model->solver()->getStrParam(OsiProbName,probName) ; |
---|
| 550 | model->solver()->activateRowCutDebugger(probName.c_str()) ; |
---|
| 551 | if (model->solver()->getRowCutDebugger()) |
---|
| 552 | std::cout << "Row cut debugger activated for " ; |
---|
| 553 | else |
---|
| 554 | std::cout << "Failed to activate row cut debugger for " ; |
---|
| 555 | std::cout << mpsName[m] << "." << std::endl ; |
---|
| 556 | } |
---|
| 557 | setCutAndHeuristicOptions(*model) ; |
---|
| 558 | /* |
---|
| 559 | More clp-specific setup. |
---|
| 560 | */ |
---|
| 561 | if (siClp) { |
---|
| 562 | # ifdef CLP_MULTIPLE_FACTORIZATIONS |
---|
| 563 | if (!modelC->factorization()->isDenseOrSmall()) { |
---|
| 564 | int denseCode = stuff ? static_cast<int> (stuff[4]) : -1; |
---|
| 565 | int smallCode = stuff ? static_cast<int> (stuff[10]) : -1; |
---|
| 566 | if (stuff && stuff[8] >= 1) { |
---|
| 567 | if (denseCode < 0) |
---|
| 568 | denseCode = 40; |
---|
| 569 | if (smallCode < 0) |
---|
| 570 | smallCode = 40; |
---|
| 571 | } |
---|
| 572 | if (denseCode > 0) |
---|
| 573 | modelC->factorization()->setGoDenseThreshold(denseCode); |
---|
| 574 | if (smallCode > 0) |
---|
| 575 | modelC->factorization()->setGoSmallThreshold(smallCode); |
---|
| 576 | if (denseCode >= modelC->numberRows()) { |
---|
| 577 | //printf("problem going dense\n"); |
---|
| 578 | //modelC->factorization()->goDenseOrSmall(modelC->numberRows()); |
---|
| 579 | } |
---|
| 580 | } |
---|
| 581 | # endif |
---|
| 582 | if (stuff && stuff[8] >= 1) { |
---|
| 583 | if (modelC->numberColumns() + modelC->numberRows() <= 10000 && |
---|
| 584 | model->fastNodeDepth() == -1) |
---|
| 585 | model->setFastNodeDepth(-9); |
---|
| 586 | } |
---|
| 587 | } |
---|
| 588 | //OsiObject * obj = new CbcBranchToFixLots(model,0.3,0.0,3,3000003); |
---|
| 589 | //model->addObjects(1,&obj); |
---|
| 590 | //delete obj; |
---|
| 591 | /* |
---|
| 592 | Finally, the actual call to solve the MIP with branch-and-cut. |
---|
| 593 | */ |
---|
| 594 | model->branchAndBound(); |
---|
| 595 | |
---|
| 596 | # ifdef CLP_FACTORIZATION_INSTRUMENT |
---|
| 597 | double facTime = factorization_instrument(0); |
---|
| 598 | std::cout << "Factorization " << facTime << " seconds." << std::endl , |
---|
| 599 | timeTakenFac += facTime; |
---|
| 600 | # endif |
---|
| 601 | |
---|
| 602 | /* |
---|
| 603 | Stage 3: Do the statistics and check the answer. |
---|
| 604 | */ |
---|
| 605 | double timeOfSolution = CoinCpuTime()+CoinCpuTimeJustChildren()-startTime; |
---|
[1286] | 606 | std::cout |
---|
[1622] | 607 | << "Cuts at root node changed objective from " |
---|
| 608 | << model->getContinuousObjective() << " to " |
---|
| 609 | << model->rootObjectiveAfterCuts() << std::endl ; |
---|
| 610 | int numberGenerators = model->numberCutGenerators(); |
---|
| 611 | for (int iGenerator = 0 ; iGenerator < numberGenerators ; iGenerator++) { |
---|
| 612 | CbcCutGenerator *generator = model->cutGenerator(iGenerator); |
---|
| 613 | # ifdef CLIQUE_ANALYSIS |
---|
| 614 | # ifndef CLP_INVESTIGATE |
---|
| 615 | CglImplication *implication = |
---|
| 616 | dynamic_cast<CglImplication*>(generator->generator()); |
---|
| 617 | if (implication) continue; |
---|
| 618 | # endif |
---|
| 619 | # endif |
---|
| 620 | std::cout |
---|
| 621 | << generator->cutGeneratorName() << " was tried " |
---|
| 622 | << generator->numberTimesEntered() << " times and created " |
---|
| 623 | << generator->numberCutsInTotal() << " cuts of which " |
---|
| 624 | << generator->numberCutsActive() |
---|
| 625 | << " were active after adding rounds of cuts"; |
---|
| 626 | if (generator->timing()) |
---|
| 627 | std::cout << " (" << generator->timeInCutGenerator() << " seconds)" ; |
---|
| 628 | std::cout << "." << std::endl; |
---|
| 629 | } |
---|
| 630 | std::cout |
---|
| 631 | << model->getNumberHeuristicSolutions() |
---|
| 632 | << " solutions found by heuristics." << std::endl ; |
---|
| 633 | int numberHeuristics = model->numberHeuristics(); |
---|
| 634 | for (int iHeuristic = 0 ; iHeuristic < numberHeuristics ; iHeuristic++) { |
---|
| 635 | CbcHeuristic *heuristic = model->heuristic(iHeuristic); |
---|
| 636 | if (heuristic->numRuns()) { |
---|
| 637 | std::cout |
---|
| 638 | << heuristic->heuristicName() << " was tried " |
---|
| 639 | << heuristic->numRuns() << " times out of " |
---|
| 640 | << heuristic->numCouldRun() << " and created " |
---|
| 641 | << heuristic->numberSolutionsFound() << " solutions." << std::endl ; |
---|
| 642 | } |
---|
| 643 | } |
---|
| 644 | /* |
---|
| 645 | Check for the correct answer. |
---|
| 646 | */ |
---|
| 647 | if (!model->status()) { |
---|
| 648 | |
---|
| 649 | double objActual = model->getObjValue() ; |
---|
| 650 | double objExpect = objValue[m] ; |
---|
| 651 | double tolerance = CoinMin(fabs(objActual),fabs(objExpect)) ; |
---|
| 652 | tolerance = CoinMax(1.0e-5,1.0e-5*tolerance) ; |
---|
| 653 | //CoinRelFltEq eq(1.0e-3) ; |
---|
| 654 | |
---|
| 655 | std::cout |
---|
| 656 | << "cbc_clp (" << mpsName[m] << ") " |
---|
| 657 | << std::setprecision(10) << objActual ; |
---|
| 658 | if (fabs(objActual-objExpect) < tolerance) { |
---|
| 659 | std::cout << std::setprecision(dfltPrecision) << "; okay" ; |
---|
| 660 | numProbSolved++; |
---|
| 661 | } else { |
---|
| 662 | std::cout |
---|
| 663 | << " != " << objExpect << std::setprecision(dfltPrecision) |
---|
| 664 | << "; error = " << fabs(objExpect-objActual) ; |
---|
| 665 | numberFailures++; |
---|
| 666 | //#ifdef COIN_DEVELOP |
---|
| 667 | //abort(); |
---|
| 668 | //#endif |
---|
| 669 | } |
---|
| 670 | } else { |
---|
| 671 | std::cout |
---|
| 672 | << "cbc_clp (" << mpsName[m] << ") status not optimal; " |
---|
| 673 | << "assuming too many nodes" ; |
---|
| 674 | } |
---|
| 675 | timeTaken += timeOfSolution; |
---|
| 676 | std::cout |
---|
| 677 | << " -- (" << model->getNodeCount() << " n / " |
---|
| 678 | << model->getIterationCount() << " i / " |
---|
[1880] | 679 | << timeOfSolution << " s) (subtotal " << timeTaken << " seconds)" |
---|
[1622] | 680 | << std::endl; |
---|
| 681 | delete model; |
---|
| 682 | } |
---|
| 683 | /* |
---|
| 684 | End main loop on MPS problems. Print a summary and calculate the return |
---|
| 685 | value. |
---|
| 686 | */ |
---|
| 687 | int returnCode = 0; |
---|
| 688 | std::cout |
---|
| 689 | << "cbc_clp solved " << numProbSolved << " out of " << numberAttempts; |
---|
| 690 | int numberOnNodes = numberAttempts-numProbSolved-numberFailures; |
---|
| 691 | if (numberFailures || numberOnNodes) { |
---|
| 692 | if (numberOnNodes) { |
---|
| 693 | std::cout << " (" << numberOnNodes << " stopped on nodes)"; |
---|
| 694 | returnCode = numberOnNodes; |
---|
| 695 | } |
---|
| 696 | if (numberFailures) { |
---|
| 697 | std::cout << " (" << numberFailures << " gave bad answer!)"; |
---|
| 698 | returnCode += 100*numberFailures; |
---|
| 699 | } |
---|
| 700 | } |
---|
| 701 | std::cout |
---|
| 702 | << " and took " << timeTaken << " seconds." << std::endl; |
---|
| 703 | |
---|
| 704 | if (testSwitch == -2) { |
---|
[1286] | 705 | if (numberFailures || numberOnNodes) { |
---|
[1622] | 706 | std::cout << "****** Unit Test failed." << std::endl ; |
---|
| 707 | std::cerr << "****** Unit Test failed." << std::endl ; |
---|
| 708 | } else { |
---|
| 709 | std::cerr << "****** Unit Test succeeded." << std::endl ; |
---|
[824] | 710 | } |
---|
[1622] | 711 | } |
---|
| 712 | # ifdef CLP_FACTORIZATION_INSTRUMENT |
---|
| 713 | std::cout |
---|
| 714 | << "Total factorization time " << timeTakenFac << "seconds." << std::endl ; |
---|
| 715 | # endif |
---|
| 716 | return (returnCode) ; |
---|
[274] | 717 | } |
---|
[1432] | 718 | |
---|