[2310] | 1 | // $Id: driverFat.cpp 1898 2013-04-09 18:06:04Z stefan $ |
---|
| 2 | // Copyright (C) 2007, International Business Machines |
---|
| 3 | // Corporation and others. All Rights Reserved. |
---|
| 4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
| 5 | |
---|
| 6 | #include <cassert> |
---|
| 7 | #include <iomanip> |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | #include "CoinPragma.hpp" |
---|
| 11 | #include "CbcModel.hpp" |
---|
| 12 | #include "CbcHeuristic.hpp" |
---|
| 13 | #include "OsiClpSolverInterface.hpp" |
---|
| 14 | #include "CbcSolver.hpp" |
---|
| 15 | #include "ClpSimplex.hpp" |
---|
| 16 | #include "ClpFactorization.hpp" |
---|
| 17 | #include "ClpPresolve.hpp" |
---|
| 18 | #include "CoinSort.hpp" |
---|
| 19 | #include "CoinHelperFunctions.hpp" |
---|
| 20 | |
---|
| 21 | #include "CoinTime.hpp" |
---|
| 22 | #include "CoinSignal.hpp" |
---|
| 23 | /* |
---|
| 24 | This shows how to trap signals. |
---|
| 25 | This just traps ctrl-c and allows user to pause and then hit S or C |
---|
| 26 | In this simple version Stop may not be effective until a heuristic has exited |
---|
| 27 | */ |
---|
| 28 | |
---|
| 29 | static CbcModel * currentBranchModel=NULL; |
---|
| 30 | extern "C" { |
---|
| 31 | static void signal_handler(int whichSignal) { |
---|
| 32 | int gotChar='X'; |
---|
| 33 | while (toupper(gotChar)!='S'&&toupper(gotChar)!='C') { |
---|
| 34 | // See what user wants to do |
---|
| 35 | fprintf(stderr,"Enter S to stop, C to continue:"); |
---|
| 36 | gotChar = getchar(); |
---|
| 37 | } |
---|
| 38 | if (currentBranchModel != NULL&&toupper(gotChar)=='S') { |
---|
| 39 | currentBranchModel->sayEventHappened(); // say why stopped |
---|
| 40 | if (currentBranchModel->heuristicModel()) |
---|
| 41 | currentBranchModel->heuristicModel()->sayEventHappened(); |
---|
| 42 | } |
---|
| 43 | return; |
---|
| 44 | } |
---|
| 45 | } |
---|
| 46 | static CoinSighandler_t saveSignal=signal(SIGINT,signal_handler); |
---|
| 47 | // Threshold below which use normal clp |
---|
| 48 | static int rowsThreshold=-1; |
---|
| 49 | //############################################################################# |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | /************************************************************************ |
---|
| 53 | |
---|
| 54 | This main program shows how to take advantage of the standalone cbc in your program, |
---|
| 55 | while still making major modifications. |
---|
| 56 | First it reads in an integer model from an mps file |
---|
| 57 | Then it initializes the integer model with cbc defaults |
---|
| 58 | Then it calls CbcMain1 passing all parameters apart from first but with callBack to modify stuff |
---|
| 59 | Finally it prints solution |
---|
| 60 | |
---|
| 61 | ************************************************************************/ |
---|
| 62 | /* Meaning of whereFrom: |
---|
| 63 | 1 after initial solve by dualsimplex etc |
---|
| 64 | 2 after preprocessing |
---|
| 65 | 3 just before branchAndBound (so user can override) |
---|
| 66 | 4 just after branchAndBound (before postprocessing) |
---|
| 67 | 5 after postprocessing |
---|
| 68 | */ |
---|
| 69 | /* Meaning of model status is as normal |
---|
| 70 | status |
---|
| 71 | -1 before branchAndBound |
---|
| 72 | 0 finished - check isProvenOptimal or isProvenInfeasible to see if solution found |
---|
| 73 | (or check value of best solution) |
---|
| 74 | 1 stopped - on maxnodes, maxsols, maxtime |
---|
| 75 | 2 difficulties so run was abandoned |
---|
| 76 | (5 event user programmed event occurred) |
---|
| 77 | |
---|
| 78 | cbc secondary status of problem |
---|
| 79 | -1 unset (status_ will also be -1) |
---|
| 80 | 0 search completed with solution |
---|
| 81 | 1 linear relaxation not feasible (or worse than cutoff) |
---|
| 82 | 2 stopped on gap |
---|
| 83 | 3 stopped on nodes |
---|
| 84 | 4 stopped on time |
---|
| 85 | 5 stopped on user event |
---|
| 86 | 6 stopped on solutions |
---|
| 87 | 7 linear relaxation unbounded |
---|
| 88 | |
---|
| 89 | but initially check if status is 0 and secondary status is 1 -> infeasible |
---|
| 90 | or you can check solver status. |
---|
| 91 | */ |
---|
| 92 | /* Return non-zero to return quickly */ |
---|
| 93 | static int callBack(CbcModel * model, int whereFrom) |
---|
| 94 | { |
---|
| 95 | int returnCode=0; |
---|
| 96 | switch (whereFrom) { |
---|
| 97 | case 1: |
---|
| 98 | case 2: |
---|
| 99 | if (!model->status()&&model->secondaryStatus()) |
---|
| 100 | returnCode=1; |
---|
| 101 | break; |
---|
| 102 | case 3: |
---|
| 103 | { |
---|
| 104 | // set up signal trapping |
---|
| 105 | saveSignal=signal(SIGINT,signal_handler); |
---|
| 106 | currentBranchModel=model; |
---|
| 107 | /******************************* |
---|
| 108 | This tells code to be normal in heuristics i.e. smaller problems |
---|
| 109 | in practice you should probably make CoinMax(...,20000) or |
---|
| 110 | some such. |
---|
| 111 | You may wish to switch off strong branching and use priorities |
---|
| 112 | or something - as strong branching uses large model |
---|
| 113 | *******************************/ |
---|
| 114 | rowsThreshold=model->getNumRows(); |
---|
| 115 | // make heuristics do more nodes |
---|
| 116 | for (int i=0;i<model->numberHeuristics();i++) { |
---|
| 117 | CbcHeuristic * heuristic = model->heuristic(i); |
---|
| 118 | heuristic->setNumberNodes(5*heuristic->numberNodes()); |
---|
| 119 | } |
---|
| 120 | // could try doing feasibility after cuts? |
---|
| 121 | //model->setSpecialOptions(33554432| |
---|
| 122 | // model->specialOptions()); |
---|
| 123 | } |
---|
| 124 | break; |
---|
| 125 | case 4: |
---|
| 126 | { |
---|
| 127 | // restore |
---|
| 128 | signal(SIGINT,saveSignal); |
---|
| 129 | currentBranchModel=NULL; |
---|
| 130 | } |
---|
| 131 | // If not good enough could skip postprocessing |
---|
| 132 | break; |
---|
| 133 | case 5: |
---|
| 134 | break; |
---|
| 135 | default: |
---|
| 136 | abort(); |
---|
| 137 | } |
---|
| 138 | return returnCode; |
---|
| 139 | } |
---|
| 140 | #include "CbcEventHandler.hpp" |
---|
| 141 | /** This is so user can trap events and do useful stuff. |
---|
| 142 | |
---|
| 143 | CbcModel model_ is available as well as anything else you care |
---|
| 144 | to pass in |
---|
| 145 | */ |
---|
| 146 | |
---|
| 147 | class MyEventHandler3 : public CbcEventHandler { |
---|
| 148 | |
---|
| 149 | public: |
---|
| 150 | /**@name Overrides */ |
---|
| 151 | //@{ |
---|
| 152 | virtual CbcAction event(CbcEvent whichEvent); |
---|
| 153 | //@} |
---|
| 154 | |
---|
| 155 | /**@name Constructors, destructor etc*/ |
---|
| 156 | //@{ |
---|
| 157 | /** Default constructor. */ |
---|
| 158 | MyEventHandler3(); |
---|
| 159 | /// Constructor with pointer to model (redundant as setEventHandler does) |
---|
| 160 | MyEventHandler3(CbcModel * model); |
---|
| 161 | /** Destructor */ |
---|
| 162 | virtual ~MyEventHandler3(); |
---|
| 163 | /** The copy constructor. */ |
---|
| 164 | MyEventHandler3(const MyEventHandler3 & rhs); |
---|
| 165 | /// Assignment |
---|
| 166 | MyEventHandler3& operator=(const MyEventHandler3 & rhs); |
---|
| 167 | /// Clone |
---|
| 168 | virtual CbcEventHandler * clone() const ; |
---|
| 169 | //@} |
---|
| 170 | |
---|
| 171 | |
---|
| 172 | protected: |
---|
| 173 | // data goes here |
---|
| 174 | }; |
---|
| 175 | //------------------------------------------------------------------- |
---|
| 176 | // Default Constructor |
---|
| 177 | //------------------------------------------------------------------- |
---|
| 178 | MyEventHandler3::MyEventHandler3 () |
---|
| 179 | : CbcEventHandler() |
---|
| 180 | { |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | //------------------------------------------------------------------- |
---|
| 184 | // Copy constructor |
---|
| 185 | //------------------------------------------------------------------- |
---|
| 186 | MyEventHandler3::MyEventHandler3 (const MyEventHandler3 & rhs) |
---|
| 187 | : CbcEventHandler(rhs) |
---|
| 188 | { |
---|
| 189 | } |
---|
| 190 | |
---|
| 191 | // Constructor with pointer to model |
---|
| 192 | MyEventHandler3::MyEventHandler3(CbcModel * model) |
---|
| 193 | : CbcEventHandler(model) |
---|
| 194 | { |
---|
| 195 | } |
---|
| 196 | |
---|
| 197 | //------------------------------------------------------------------- |
---|
| 198 | // Destructor |
---|
| 199 | //------------------------------------------------------------------- |
---|
| 200 | MyEventHandler3::~MyEventHandler3 () |
---|
| 201 | { |
---|
| 202 | } |
---|
| 203 | |
---|
| 204 | //---------------------------------------------------------------- |
---|
| 205 | // Assignment operator |
---|
| 206 | //------------------------------------------------------------------- |
---|
| 207 | MyEventHandler3 & |
---|
| 208 | MyEventHandler3::operator=(const MyEventHandler3& rhs) |
---|
| 209 | { |
---|
| 210 | if (this != &rhs) { |
---|
| 211 | CbcEventHandler::operator=(rhs); |
---|
| 212 | } |
---|
| 213 | return *this; |
---|
| 214 | } |
---|
| 215 | //------------------------------------------------------------------- |
---|
| 216 | // Clone |
---|
| 217 | //------------------------------------------------------------------- |
---|
| 218 | CbcEventHandler * MyEventHandler3::clone() const |
---|
| 219 | { |
---|
| 220 | return new MyEventHandler3(*this); |
---|
| 221 | } |
---|
| 222 | |
---|
| 223 | CbcEventHandler::CbcAction |
---|
| 224 | MyEventHandler3::event(CbcEvent whichEvent) |
---|
| 225 | { |
---|
| 226 | // If in sub tree carry on |
---|
| 227 | if (!model_->parentModel()) { |
---|
| 228 | if (whichEvent==solution||whichEvent==heuristicSolution) { |
---|
| 229 | #ifdef STOP_EARLY |
---|
| 230 | return stop; // say finished |
---|
| 231 | #else |
---|
| 232 | // If preprocessing was done solution will be to processed model |
---|
| 233 | #if 0 |
---|
| 234 | int numberColumns = model_->getNumCols(); |
---|
| 235 | const double * bestSolution = model_->bestSolution(); |
---|
| 236 | assert (bestSolution); |
---|
| 237 | printf("value of solution is %g\n",model_->getObjValue()); |
---|
| 238 | for (int i=0;i<numberColumns;i++) { |
---|
| 239 | if (fabs(bestSolution[i])>1.0e-8) |
---|
| 240 | printf("%d %g\n",i,bestSolution[i]); |
---|
| 241 | } |
---|
| 242 | #endif |
---|
| 243 | return noAction; // carry on |
---|
| 244 | #endif |
---|
| 245 | } else { |
---|
| 246 | return noAction; // carry on |
---|
| 247 | } |
---|
| 248 | } else { |
---|
| 249 | return noAction; // carry on |
---|
| 250 | } |
---|
| 251 | } |
---|
| 252 | //############################################################################# |
---|
| 253 | |
---|
| 254 | /** |
---|
| 255 | |
---|
| 256 | This is to allow the user to replace initialSolve and resolve |
---|
| 257 | */ |
---|
| 258 | |
---|
| 259 | class CbcSolverShortFat : public OsiClpSolverInterface { |
---|
| 260 | |
---|
| 261 | public: |
---|
| 262 | //--------------------------------------------------------------------------- |
---|
| 263 | /**@name Solve methods */ |
---|
| 264 | //@{ |
---|
| 265 | /// Solve initial LP relaxation |
---|
| 266 | virtual void initialSolve(); |
---|
| 267 | |
---|
| 268 | /// Resolve an LP relaxation after problem modification |
---|
| 269 | virtual void resolve(); |
---|
| 270 | |
---|
| 271 | //@} |
---|
| 272 | |
---|
| 273 | |
---|
| 274 | /**@name Constructors and destructors */ |
---|
| 275 | //@{ |
---|
| 276 | /// Default Constructor |
---|
| 277 | CbcSolverShortFat (); |
---|
| 278 | |
---|
| 279 | /// Clone |
---|
| 280 | virtual OsiSolverInterface * clone(bool CopyData=true) const; |
---|
| 281 | |
---|
| 282 | /// Copy constructor |
---|
| 283 | CbcSolverShortFat (const CbcSolverShortFat &); |
---|
| 284 | |
---|
| 285 | /// Assignment operator |
---|
| 286 | CbcSolverShortFat & operator=(const CbcSolverShortFat& rhs); |
---|
| 287 | |
---|
| 288 | /// Destructor |
---|
| 289 | virtual ~CbcSolverShortFat (); |
---|
| 290 | |
---|
| 291 | //@} |
---|
| 292 | |
---|
| 293 | |
---|
| 294 | |
---|
| 295 | //--------------------------------------------------------------------------- |
---|
| 296 | |
---|
| 297 | private: |
---|
| 298 | |
---|
| 299 | /**@name Private member data */ |
---|
| 300 | //@{ |
---|
| 301 | //@} |
---|
| 302 | }; |
---|
| 303 | static bool firstSolve=true; |
---|
| 304 | void CbcSolverShortFat::initialSolve() |
---|
| 305 | { |
---|
| 306 | ClpSimplex * model = getModelPtr(); |
---|
| 307 | if (model->numberRows()<rowsThreshold) { |
---|
| 308 | OsiClpSolverInterface::initialSolve(); |
---|
| 309 | return; |
---|
| 310 | } |
---|
| 311 | #if LOGLEVEL>1 |
---|
| 312 | double time1=CoinCpuTime(); |
---|
| 313 | #endif |
---|
| 314 | ClpPresolve pinfo; |
---|
| 315 | #define INCREASE 6 |
---|
| 316 | #define PERTURB 5 |
---|
| 317 | #ifdef PERTURB |
---|
| 318 | bool externalPerturb=true; |
---|
| 319 | #endif |
---|
| 320 | if (firstSolve) { // before preprocessing |
---|
| 321 | model = pinfo.presolvedModel(*model, 1.0e-8, false, 5, false); |
---|
| 322 | } else { |
---|
| 323 | externalPerturb=false; |
---|
| 324 | /* do initial factorization to get infeasibilities |
---|
| 325 | maybe fix all non basic |
---|
| 326 | initial fix checks if Osi already has basis - yes it has |
---|
| 327 | |
---|
| 328 | */ |
---|
| 329 | setBasis(basis_,model); |
---|
| 330 | } |
---|
| 331 | #define LOGLEVEL 0 |
---|
| 332 | #if LOGLEVEL<3 |
---|
| 333 | model->setLogLevel(0); |
---|
| 334 | #endif |
---|
| 335 | int numberColumns = model->numberColumns(); |
---|
| 336 | int originalNumberRows = model->numberRows(); |
---|
| 337 | // change factorization frequency from 200 |
---|
| 338 | model->setFactorizationFrequency(100 + model->numberRows() / 50); |
---|
| 339 | int numberIterations=0; |
---|
| 340 | #ifdef PERTURB |
---|
| 341 | double * objective = model->objective(); |
---|
| 342 | double * saveObjective = NULL; |
---|
| 343 | if (externalPerturb) { |
---|
| 344 | saveObjective = CoinCopyOfArray(objective,numberColumns); |
---|
| 345 | double multiplier=1.0; |
---|
| 346 | for (int i=0;i<PERTURB;i++) |
---|
| 347 | multiplier *= 0.1; |
---|
| 348 | for (int i=0;i<numberColumns;i++) { |
---|
| 349 | double value = CoinDrand48()*multiplier; |
---|
| 350 | // should be more sophisticated |
---|
| 351 | if (value>1.0e-7) { |
---|
| 352 | if (objective[i]<0.0) |
---|
| 353 | value = -value; |
---|
| 354 | objective[i] += value; |
---|
| 355 | } |
---|
| 356 | } |
---|
| 357 | } |
---|
| 358 | #endif |
---|
| 359 | |
---|
| 360 | // We will need arrays to choose rows to add |
---|
| 361 | double * weight = new double [originalNumberRows]; |
---|
| 362 | int * sort = new int [originalNumberRows]; |
---|
| 363 | int numberSort = 0; |
---|
| 364 | char * take = new char [originalNumberRows]; |
---|
| 365 | |
---|
| 366 | const double * rowLower = model->rowLower(); |
---|
| 367 | const double * rowUpper = model->rowUpper(); |
---|
| 368 | int iRow, iColumn; |
---|
| 369 | // Set up initial list |
---|
| 370 | numberSort = 0; |
---|
| 371 | int numberNonBasicRows=0; |
---|
| 372 | for (iRow = 0; iRow < originalNumberRows; iRow++) { |
---|
| 373 | weight[iRow] = 1.123e50; |
---|
| 374 | if (model->getRowStatus(iRow)!=ClpSimplex::basic) { |
---|
| 375 | sort[numberSort++] = iRow; |
---|
| 376 | weight[iRow] = -10.0; |
---|
| 377 | numberNonBasicRows++; |
---|
| 378 | } else if (rowLower[iRow] == rowUpper[iRow]) { |
---|
| 379 | sort[numberSort++] = iRow; |
---|
| 380 | weight[iRow] = 0.0; |
---|
| 381 | } |
---|
| 382 | } |
---|
| 383 | numberSort /= 2; |
---|
| 384 | numberSort = CoinMax(numberSort,numberNonBasicRows); |
---|
| 385 | // Just add this number of rows each time in small problem |
---|
| 386 | int smallNumberRows = 2 * numberColumns; |
---|
| 387 | smallNumberRows = CoinMin(smallNumberRows, originalNumberRows / 20); |
---|
| 388 | // and pad out with random rows |
---|
| 389 | double ratio = (static_cast<double>(smallNumberRows - numberSort)) / |
---|
| 390 | (static_cast<double>(originalNumberRows)); |
---|
| 391 | bool primalInfeasible=false; |
---|
| 392 | if (firstSolve) { |
---|
| 393 | for (iRow = 0; iRow < originalNumberRows; iRow++) { |
---|
| 394 | if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio) |
---|
| 395 | sort[numberSort++] = iRow; |
---|
| 396 | } |
---|
| 397 | } |
---|
| 398 | /* This is optional. |
---|
| 399 | The best thing to do is to miss out random rows and do a set which makes dual feasible. |
---|
| 400 | If that is not possible then make sure variables have bounds. |
---|
| 401 | |
---|
| 402 | One way that normally works is to automatically tighten bounds. |
---|
| 403 | */ |
---|
| 404 | if (firstSolve) { |
---|
| 405 | // However for some we need to do anyway |
---|
| 406 | double * columnLower = model->columnLower(); |
---|
| 407 | double * columnUpper = model->columnUpper(); |
---|
| 408 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 409 | columnLower[iColumn] = CoinMax(-1.0e6, columnLower[iColumn]); |
---|
| 410 | columnUpper[iColumn] = CoinMin(1.0e6, columnUpper[iColumn]); |
---|
| 411 | } |
---|
| 412 | } |
---|
| 413 | double * fullSolution = model->primalRowSolution(); |
---|
| 414 | |
---|
| 415 | // Just do this number of passes |
---|
| 416 | int maxPass = 50; |
---|
| 417 | // And take out slack rows until this pass |
---|
| 418 | int takeOutPass = INCREASE; |
---|
| 419 | int iPass; |
---|
| 420 | |
---|
[2346] | 421 | const CoinBigIndex * start = model->clpMatrix()->getVectorStarts(); |
---|
[2310] | 422 | const int * length = model->clpMatrix()->getVectorLengths(); |
---|
| 423 | const int * row = model->clpMatrix()->getIndices(); |
---|
| 424 | int * whichColumns = new int [numberColumns]; |
---|
| 425 | for (int iRow = 0; iRow < originalNumberRows; iRow++) |
---|
| 426 | weight[iRow]=0.0; |
---|
| 427 | for (int i = 0; i < numberSort; i++) { |
---|
| 428 | int iRow=sort[i]; |
---|
| 429 | weight[iRow]=1.0; |
---|
| 430 | } |
---|
| 431 | |
---|
| 432 | int numberSmallColumns = 0; |
---|
| 433 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 434 | bool take=false;; |
---|
| 435 | for (int j=start[iColumn];j<start[iColumn]+length[iColumn];j++) { |
---|
| 436 | int iRow=row[j]; |
---|
| 437 | if (weight[iRow]) { |
---|
| 438 | take=true; |
---|
| 439 | break; |
---|
| 440 | } |
---|
| 441 | } |
---|
| 442 | if (take) |
---|
| 443 | whichColumns[numberSmallColumns++] = iColumn; |
---|
| 444 | } |
---|
| 445 | #if LOGLEVEL>1 |
---|
| 446 | printf("%d rows, %d columns in initial problem\n", numberSort,numberSmallColumns); |
---|
| 447 | #endif |
---|
| 448 | for (iPass = 0; iPass < maxPass; iPass++) { |
---|
| 449 | #if LOGLEVEL>2 |
---|
| 450 | printf("Start of pass %d\n", iPass); |
---|
| 451 | #endif |
---|
| 452 | // Cleaner this way |
---|
| 453 | std::sort(sort, sort + numberSort); |
---|
| 454 | // Create small problem |
---|
| 455 | ClpSimplex small(model, numberSort, sort, numberSmallColumns, whichColumns); |
---|
| 456 | small.setFactorizationFrequency(100 + numberSort / 200); |
---|
| 457 | #ifndef PERTURB |
---|
| 458 | small.setPerturbation(50); |
---|
| 459 | #else |
---|
| 460 | if (!externalPerturb) |
---|
| 461 | small.setPerturbation(50); |
---|
| 462 | #endif |
---|
| 463 | #if LOGLEVEL>2 |
---|
| 464 | small.setLogLevel(1); |
---|
| 465 | #endif |
---|
| 466 | // A variation is to just do N iterations |
---|
| 467 | //if (iPass) |
---|
| 468 | //small.setMaximumIterations(100); |
---|
| 469 | // Solve |
---|
| 470 | //small.factorization()->messageLevel(8); |
---|
| 471 | small.dual(); |
---|
| 472 | numberIterations+= small.numberIterations(); |
---|
| 473 | primalInfeasible = (small.status() ==1); |
---|
| 474 | if (primalInfeasible) |
---|
| 475 | break; |
---|
| 476 | bool dualInfeasible = (small.status() == 2); |
---|
| 477 | // move solution back |
---|
| 478 | double * solution = model->primalColumnSolution(); |
---|
| 479 | const double * smallSolution = small.primalColumnSolution(); |
---|
| 480 | for (int j = 0; j < numberSmallColumns; j++) { |
---|
| 481 | iColumn = whichColumns[j]; |
---|
| 482 | solution[iColumn] = smallSolution[j]; |
---|
| 483 | model->setColumnStatus(iColumn, small.getColumnStatus(j)); |
---|
| 484 | } |
---|
| 485 | for (iRow = 0; iRow < numberSort; iRow++) { |
---|
| 486 | int kRow = sort[iRow]; |
---|
| 487 | model->setRowStatus(kRow, small.getRowStatus(iRow)); |
---|
| 488 | } |
---|
| 489 | // compute full solution |
---|
| 490 | memset(fullSolution, 0, originalNumberRows * sizeof(double)); |
---|
| 491 | model->clpMatrix()->times(1.0, model->primalColumnSolution(), fullSolution); |
---|
| 492 | if (iPass != maxPass - 1) { |
---|
| 493 | // Mark row as not looked at |
---|
| 494 | for (iRow = 0; iRow < originalNumberRows; iRow++) |
---|
| 495 | weight[iRow] = 1.123e50; |
---|
| 496 | // Look at rows already in small problem |
---|
| 497 | int iSort; |
---|
| 498 | int numberDropped = 0; |
---|
| 499 | int numberKept = 0; |
---|
| 500 | int numberBinding = 0; |
---|
| 501 | int numberInfeasibilities = 0; |
---|
| 502 | double sumInfeasibilities = 0.0; |
---|
| 503 | for (iSort = 0; iSort < numberSort; iSort++) { |
---|
| 504 | iRow = sort[iSort]; |
---|
| 505 | if (model->getRowStatus(iRow) == ClpSimplex::basic) { |
---|
| 506 | // Basic - we can get rid of if early on |
---|
| 507 | if (iPass < takeOutPass && !dualInfeasible) { |
---|
| 508 | // may have hit max iterations so check |
---|
| 509 | double infeasibility = CoinMax(fullSolution[iRow] - rowUpper[iRow], |
---|
| 510 | rowLower[iRow] - fullSolution[iRow]); |
---|
| 511 | weight[iRow] = -infeasibility; |
---|
| 512 | if (infeasibility > 1.0e-8) { |
---|
| 513 | numberInfeasibilities++; |
---|
| 514 | sumInfeasibilities += infeasibility; |
---|
| 515 | } else { |
---|
| 516 | weight[iRow] = 1.0; |
---|
| 517 | numberDropped++; |
---|
| 518 | } |
---|
| 519 | } else { |
---|
| 520 | // keep |
---|
| 521 | weight[iRow] = -1.0e40; |
---|
| 522 | numberKept++; |
---|
| 523 | } |
---|
| 524 | } else { |
---|
| 525 | // keep |
---|
| 526 | weight[iRow] = -1.0e50; |
---|
| 527 | numberKept++; |
---|
| 528 | numberBinding++; |
---|
| 529 | } |
---|
| 530 | } |
---|
| 531 | // Now rest |
---|
| 532 | for (iRow = 0; iRow < originalNumberRows; iRow++) { |
---|
| 533 | sort[iRow] = iRow; |
---|
| 534 | if (weight[iRow] == 1.123e50) { |
---|
| 535 | // not looked at yet |
---|
| 536 | double infeasibility = CoinMax(fullSolution[iRow] - rowUpper[iRow], |
---|
| 537 | rowLower[iRow] - fullSolution[iRow]); |
---|
| 538 | weight[iRow] = -infeasibility; |
---|
| 539 | if (infeasibility > 1.0e-8) { |
---|
| 540 | numberInfeasibilities++; |
---|
| 541 | sumInfeasibilities += infeasibility; |
---|
| 542 | } |
---|
| 543 | } |
---|
| 544 | } |
---|
| 545 | // sort |
---|
| 546 | CoinSort_2(weight, weight + originalNumberRows, sort); |
---|
| 547 | numberSort = CoinMin(originalNumberRows, smallNumberRows + numberKept); |
---|
| 548 | memset(take, 0, originalNumberRows); |
---|
| 549 | for (iRow = 0; iRow < numberSort; iRow++) |
---|
| 550 | take[sort[iRow]] = 1; |
---|
| 551 | numberSmallColumns = 0; |
---|
| 552 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 553 | int n = 0; |
---|
| 554 | for (int j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) { |
---|
| 555 | int iRow = row[j]; |
---|
| 556 | if (take[iRow]) |
---|
| 557 | n++; |
---|
| 558 | } |
---|
| 559 | if (n) |
---|
| 560 | whichColumns[numberSmallColumns++] = iColumn; |
---|
| 561 | } |
---|
| 562 | #if LOGLEVEL>1 |
---|
| 563 | printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n", |
---|
| 564 | numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns); |
---|
| 565 | printf("%d rows are infeasible - sum is %g\n", numberInfeasibilities, |
---|
| 566 | sumInfeasibilities); |
---|
| 567 | #endif |
---|
| 568 | if (!numberInfeasibilities) { |
---|
| 569 | #if LOGLEVEL>1 |
---|
| 570 | printf("Exiting as looks optimal\n"); |
---|
| 571 | #endif |
---|
| 572 | break; |
---|
| 573 | } |
---|
| 574 | numberInfeasibilities = 0; |
---|
| 575 | sumInfeasibilities = 0.0; |
---|
| 576 | for (iSort = 0; iSort < numberSort; iSort++) { |
---|
| 577 | if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) { |
---|
| 578 | numberInfeasibilities++; |
---|
| 579 | sumInfeasibilities += -weight[iSort]; |
---|
| 580 | } |
---|
| 581 | } |
---|
| 582 | #if LOGLEVEL>1 |
---|
| 583 | printf("in small model %d rows are infeasible - sum is %g\n", numberInfeasibilities, |
---|
| 584 | sumInfeasibilities); |
---|
| 585 | #endif |
---|
| 586 | } |
---|
| 587 | } |
---|
| 588 | delete [] weight; |
---|
| 589 | delete [] sort; |
---|
| 590 | delete [] whichColumns; |
---|
| 591 | delete [] take; |
---|
| 592 | #ifdef PERTURB |
---|
| 593 | if (externalPerturb) { |
---|
| 594 | memcpy(objective,saveObjective,numberColumns*sizeof(double)); |
---|
| 595 | delete [] saveObjective; |
---|
| 596 | } |
---|
| 597 | model->setPerturbation(50); |
---|
| 598 | #endif |
---|
| 599 | if(!primalInfeasible) { |
---|
| 600 | model->primal(1); |
---|
| 601 | numberIterations+= model->numberIterations(); |
---|
| 602 | } else { |
---|
| 603 | model->setProblemStatus(1); |
---|
| 604 | } |
---|
| 605 | model->setNumberIterations(numberIterations); |
---|
| 606 | if (firstSolve) { |
---|
| 607 | pinfo.postsolve(true); |
---|
| 608 | model=getModelPtr(); |
---|
| 609 | model->primal(1); |
---|
| 610 | firstSolve=false; |
---|
| 611 | } |
---|
| 612 | basis_ = getBasis(model); |
---|
| 613 | #if LOGLEVEL>1 |
---|
| 614 | printf("solve took %g seconds and %d iterations\n",CoinCpuTime()-time1, |
---|
| 615 | numberIterations); |
---|
| 616 | #endif |
---|
| 617 | } |
---|
| 618 | |
---|
| 619 | //----------------------------------------------------------------------------- |
---|
| 620 | void CbcSolverShortFat::resolve() |
---|
| 621 | { |
---|
| 622 | ClpSimplex * model = getModelPtr(); |
---|
| 623 | if (model->numberRows()<rowsThreshold) { |
---|
| 624 | OsiClpSolverInterface::resolve(); |
---|
| 625 | } else { |
---|
| 626 | #if LOGLEVEL>1 |
---|
| 627 | printf("resolve\n"); |
---|
| 628 | #endif |
---|
| 629 | initialSolve(); |
---|
| 630 | } |
---|
| 631 | return; |
---|
| 632 | } |
---|
| 633 | |
---|
| 634 | //############################################################################# |
---|
| 635 | // Constructors, destructors clone and assignment |
---|
| 636 | //############################################################################# |
---|
| 637 | |
---|
| 638 | //------------------------------------------------------------------- |
---|
| 639 | // Default Constructor |
---|
| 640 | //------------------------------------------------------------------- |
---|
| 641 | CbcSolverShortFat::CbcSolverShortFat () |
---|
| 642 | : OsiClpSolverInterface() |
---|
| 643 | { |
---|
| 644 | } |
---|
| 645 | |
---|
| 646 | //------------------------------------------------------------------- |
---|
| 647 | // Clone |
---|
| 648 | //------------------------------------------------------------------- |
---|
| 649 | OsiSolverInterface * |
---|
| 650 | CbcSolverShortFat::clone(bool CopyData) const |
---|
| 651 | { |
---|
| 652 | if (CopyData) { |
---|
| 653 | return new CbcSolverShortFat(*this); |
---|
| 654 | } else { |
---|
| 655 | printf("warning CbcSolveUser clone with copyData false\n"); |
---|
| 656 | return new CbcSolverShortFat(); |
---|
| 657 | } |
---|
| 658 | } |
---|
| 659 | |
---|
| 660 | |
---|
| 661 | //------------------------------------------------------------------- |
---|
| 662 | // Copy constructor |
---|
| 663 | //------------------------------------------------------------------- |
---|
| 664 | CbcSolverShortFat::CbcSolverShortFat ( |
---|
| 665 | const CbcSolverShortFat & rhs) |
---|
| 666 | : OsiClpSolverInterface(rhs) |
---|
| 667 | { |
---|
| 668 | } |
---|
| 669 | |
---|
| 670 | //------------------------------------------------------------------- |
---|
| 671 | // Destructor |
---|
| 672 | //------------------------------------------------------------------- |
---|
| 673 | CbcSolverShortFat::~CbcSolverShortFat () |
---|
| 674 | { |
---|
| 675 | } |
---|
| 676 | |
---|
| 677 | //------------------------------------------------------------------- |
---|
| 678 | // Assignment operator |
---|
| 679 | //------------------------------------------------------------------- |
---|
| 680 | CbcSolverShortFat & |
---|
| 681 | CbcSolverShortFat::operator=(const CbcSolverShortFat& rhs) |
---|
| 682 | { |
---|
| 683 | if (this != &rhs) { |
---|
| 684 | OsiClpSolverInterface::operator=(rhs); |
---|
| 685 | } |
---|
| 686 | return *this; |
---|
| 687 | } |
---|
| 688 | //------------------------------------------------------------------- |
---|
| 689 | |
---|
| 690 | int main (int argc, const char *argv[]) |
---|
| 691 | { |
---|
| 692 | |
---|
| 693 | CbcSolverShortFat solver1; |
---|
| 694 | // Read in model using argv[1] |
---|
| 695 | // and assert that it is a clean model |
---|
| 696 | std::string mpsFileName; |
---|
| 697 | if (argc < 2) { |
---|
| 698 | fprintf(stderr, "Do not know where to find input file.\n"); |
---|
| 699 | exit(1); |
---|
| 700 | } |
---|
| 701 | if (argc>=2) mpsFileName = argv[1]; |
---|
| 702 | int numMpsReadErrors; |
---|
| 703 | if (!strstr(mpsFileName.c_str(),".lp")) |
---|
| 704 | numMpsReadErrors = solver1.readMps(mpsFileName.c_str(),""); |
---|
| 705 | else |
---|
| 706 | numMpsReadErrors = solver1.readLp(mpsFileName.c_str(),1.0e-12); |
---|
| 707 | if( numMpsReadErrors != 0 ) |
---|
| 708 | { |
---|
| 709 | printf("%d errors reading MPS file\n", numMpsReadErrors); |
---|
| 710 | return numMpsReadErrors; |
---|
| 711 | } |
---|
| 712 | // Tell solver to return fast if presolve or initial solve infeasible |
---|
| 713 | solver1.getModelPtr()->setMoreSpecialOptions(3); |
---|
| 714 | |
---|
| 715 | // Pass to Cbc initialize defaults |
---|
| 716 | CbcModel modelA(solver1); |
---|
| 717 | CbcModel * model = &modelA; |
---|
| 718 | CbcMain0(modelA); |
---|
| 719 | // Event handler |
---|
| 720 | MyEventHandler3 eventHandler; |
---|
| 721 | model->passInEventHandler(&eventHandler); |
---|
| 722 | /* Now go into code for standalone solver |
---|
| 723 | Could copy arguments and add -quit at end to be safe |
---|
| 724 | but this will do |
---|
| 725 | */ |
---|
| 726 | if (argc>2) { |
---|
| 727 | CbcMain1(argc-1,argv+1,modelA,callBack); |
---|
| 728 | } else { |
---|
| 729 | const char * argv2[]={"driverFat","-solve","-quit"}; |
---|
| 730 | CbcMain1(3,argv2,modelA,callBack); |
---|
| 731 | } |
---|
| 732 | // Solver was cloned so get current copy |
---|
| 733 | OsiSolverInterface * solver = model->solver(); |
---|
| 734 | // Print solution if finished (could get from model->bestSolution() as well |
---|
| 735 | |
---|
| 736 | if (model->bestSolution()) { |
---|
| 737 | |
---|
| 738 | const double * solution = solver->getColSolution(); |
---|
| 739 | |
---|
| 740 | int iColumn; |
---|
| 741 | int numberColumns = solver->getNumCols(); |
---|
| 742 | std::cout<<std::setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14); |
---|
| 743 | |
---|
| 744 | std::cout<<"--------------------------------------"<<std::endl; |
---|
| 745 | // names may not be in current solver - use original |
---|
| 746 | |
---|
| 747 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
| 748 | double value=solution[iColumn]; |
---|
| 749 | if (fabs(value)>1.0e-7&&solver->isInteger(iColumn)) |
---|
| 750 | std::cout<<std::setw(6)<<iColumn<<" "<<std::setw(8)<<setiosflags(std::ios::left)<<solver1.getModelPtr()->columnName(iColumn) |
---|
| 751 | <<resetiosflags(std::ios::adjustfield)<<std::setw(14)<<" "<<value<<std::endl; |
---|
| 752 | } |
---|
| 753 | std::cout<<"--------------------------------------"<<std::endl; |
---|
| 754 | |
---|
| 755 | std::cout<<std::resetiosflags(std::ios::fixed|std::ios::showpoint|std::ios::scientific); |
---|
| 756 | } else { |
---|
| 757 | std::cout<<" No solution!"<<std::endl; |
---|
| 758 | } |
---|
| 759 | return 0; |
---|
| 760 | } |
---|