Changeset 2469 for trunk/Cbc/examples/CbcSolver2.cpp
 Timestamp:
 Jan 6, 2019 6:17:46 PM (9 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Cbc/examples/CbcSolver2.cpp
r1574 r2469 18 18 #include "CoinModel.hpp" 19 19 20 static int timesBad_ =0;21 static int iterationsBad_ =0;20 static int timesBad_ = 0; 21 static int iterationsBad_ = 0; 22 22 //############################################################################# 23 23 // Solve methods … … 26 26 { 27 27 modelPtr_>scaling(0); 28 setBasis(basis_, modelPtr_);28 setBasis(basis_, modelPtr_); 29 29 // Do long thin by sprint 30 30 ClpSolve options; 31 31 options.setSolveType(ClpSolve::usePrimalorSprint); 32 32 options.setPresolveType(ClpSolve::presolveOff); 33 options.setSpecialOption(1, 3,30);33 options.setSpecialOption(1, 3, 30); 34 34 modelPtr_>initialSolve(options); 35 35 basis_ = getBasis(modelPtr_); … … 41 41 { 42 42 int numberColumns = modelPtr_>numberColumns(); 43 if ((count_ <10&&algorithm_==2)!algorithm_) {43 if ((count_ < 10 && algorithm_ == 2)  !algorithm_) { 44 44 OsiClpSolverInterface::resolve(); 45 if (modelPtr_>status() ==0) {45 if (modelPtr_>status() == 0) { 46 46 count_++; 47 double * 47 double *solution = modelPtr_>primalColumnSolution(); 48 48 int i; 49 for (i =0;i<numberColumns;i++) {50 if (solution[i]>1.0e6modelPtr_>getStatus(i)==ClpSimplex::basic) {51 node_[i]=CoinMax(count_,node_[i]);52 53 49 for (i = 0; i < numberColumns; i++) { 50 if (solution[i] > 1.0e6  modelPtr_>getStatus(i) == ClpSimplex::basic) { 51 node_[i] = CoinMax(count_, node_[i]); 52 howMany_[i]++; 53 } 54 54 } 55 55 } else { 56 if (!algorithm_ ==2)57 56 if (!algorithm_ == 2) 57 printf("infeasible early on\n"); 58 58 } 59 59 } else { 60 60 // use counts 61 int numberRows =modelPtr_>numberRows();62 int * 63 int * whichColumn = new int[numberColumns];61 int numberRows = modelPtr_>numberRows(); 62 int *whichRow = new int[numberRows]; 63 int *whichColumn = new int[numberColumns]; 64 64 int i; 65 const double * 66 const double * 67 const double * 68 bool equality =false;69 for (i =0;i<numberRows;i++) {70 if (rowUpper[i] ==1.0) {71 equality =true;65 const double *lower = modelPtr_>columnLower(); 66 const double *upper = modelPtr_>columnUpper(); 67 const double *rowUpper = modelPtr_>rowUpper(); 68 bool equality = false; 69 for (i = 0; i < numberRows; i++) { 70 if (rowUpper[i] == 1.0) { 71 equality = true; 72 72 break; 73 73 } 74 74 } 75 setBasis(basis_, modelPtr_);76 int nNewCol =0;75 setBasis(basis_, modelPtr_); 76 int nNewCol = 0; 77 77 // Column copy 78 78 //const double * element = modelPtr_>matrix()>getElements(); 79 const int * 80 const CoinBigIndex * 81 const int * 82 83 int * 84 memset(rowActivity, 0,numberRows*sizeof(int));85 int * 86 memset(rowActivity2, 0,numberRows*sizeof(int));87 char * 88 memset(mark, 0,numberColumns);79 const int *row = modelPtr_>matrix()>getIndices(); 80 const CoinBigIndex *columnStart = modelPtr_>matrix()>getVectorStarts(); 81 const int *columnLength = modelPtr_>matrix()>getVectorLengths(); 82 83 int *rowActivity = new int[numberRows]; 84 memset(rowActivity, 0, numberRows * sizeof(int)); 85 int *rowActivity2 = new int[numberRows]; 86 memset(rowActivity2, 0, numberRows * sizeof(int)); 87 char *mark = new char[numberColumns]; 88 memset(mark, 0, numberColumns); 89 89 // Get rows which are satisfied 90 for (i =0;i<numberColumns;i++) {91 if (lower[i] >0.0) {90 for (i = 0; i < numberColumns; i++) { 91 if (lower[i] > 0.0) { 92 92 CoinBigIndex j; 93 for (j =columnStart[i];94 j <columnStart[i]+columnLength[i];j++) {95 int iRow =row[j];96 rowActivity2[iRow] 93 for (j = columnStart[i]; 94 j < columnStart[i] + columnLength[i]; j++) { 95 int iRow = row[j]; 96 rowActivity2[iRow]++; 97 97 } 98 98 } else if (!upper[i]) { 99 mark[i] =2; // no good99 mark[i] = 2; // no good 100 100 } 101 101 } 102 102 // If equality  check not infeasible 103 103 if (equality) { 104 bool feasible =true;105 for (i =0;i<numberRows;i++) {106 if (rowActivity2[i] >1) {107 feasible =false;104 bool feasible = true; 105 for (i = 0; i < numberRows; i++) { 106 if (rowActivity2[i] > 1) { 107 feasible = false; 108 108 break; 109 109 } 110 110 } 111 111 if (!feasible) { 112 delete 113 delete 112 delete[] rowActivity; 113 delete[] rowActivity2; 114 114 modelPtr_>setProblemStatus(1); 115 delete 116 delete 117 delete 115 delete[] whichRow; 116 delete[] whichColumn; 117 delete[] mark; 118 118 printf("infeasible by inspection (over)\n"); 119 119 return; 120 120 } 121 121 } 122 int nNoGood =0;123 for (i =0;i<numberColumns;i++) {124 if (mark[i] ==2) {122 int nNoGood = 0; 123 for (i = 0; i < numberColumns; i++) { 124 if (mark[i] == 2) { 125 125 nNoGood++; 126 126 continue; 127 127 } 128 128 bool choose; 129 if (algorithm_ ==1)129 if (algorithm_ == 1) 130 130 choose = true; 131 131 else 132 choose = (node_[i] >count_memory_&&node_[i]>0);132 choose = (node_[i] > count_  memory_ && node_[i] > 0); 133 133 bool any; 134 134 if (equality) { 135 135 // See if forced to be zero 136 136 CoinBigIndex j; 137 any =true;138 for (j =columnStart[i];139 j <columnStart[i]+columnLength[i];j++) {140 int iRow =row[j];137 any = true; 138 for (j = columnStart[i]; 139 j < columnStart[i] + columnLength[i]; j++) { 140 int iRow = row[j]; 141 141 if (rowActivity2[iRow]) 142 any =false; // can't be in142 any = false; // can't be in 143 143 } 144 144 } else { 145 145 // See if not useful 146 146 CoinBigIndex j; 147 any =false;148 for (j =columnStart[i];149 j <columnStart[i]+columnLength[i];j++) {150 int iRow =row[j];147 any = false; 148 for (j = columnStart[i]; 149 j < columnStart[i] + columnLength[i]; j++) { 150 int iRow = row[j]; 151 151 if (!rowActivity2[iRow]) 152 any =true; // useful153 } 154 } 155 if (!any &&!lower[i]) {156 choose =false;152 any = true; // useful 153 } 154 } 155 if (!any && !lower[i]) { 156 choose = false; 157 157 // and say can't be useful 158 mark[i] =2;158 mark[i] = 2; 159 159 nNoGood++; 160 160 } 161 if (strategy_ &&modelPtr_>getColumnStatus(i)==ClpSimplex::basic)162 choose =true;163 if (choose lower[i]>0.0) {164 mark[i] =1;165 whichColumn[nNewCol++]=i;161 if (strategy_ && modelPtr_>getColumnStatus(i) == ClpSimplex::basic) 162 choose = true; 163 if (choose  lower[i] > 0.0) { 164 mark[i] = 1; 165 whichColumn[nNewCol++] = i; 166 166 CoinBigIndex j; 167 167 double value = upper[i]; 168 168 if (value) { 169 for (j =columnStart[i];170 j <columnStart[i]+columnLength[i];j++) {171 int iRow =row[j];172 rowActivity[iRow] 169 for (j = columnStart[i]; 170 j < columnStart[i] + columnLength[i]; j++) { 171 int iRow = row[j]; 172 rowActivity[iRow]++; 173 173 } 174 174 } … … 178 178 CoinModel build; 179 179 if (equality) { 180 int row =0;181 for (i =0;i<numberRows;i++) {180 int row = 0; 181 for (i = 0; i < numberRows; i++) { 182 182 // put in all rows if wanted 183 if (strategy_)184 rowActivity2[i] =0;183 if (strategy_) 184 rowActivity2[i] = 0; 185 185 if (!rowActivity2[i]) { 186 double element =1.0;187 build.addColumn(1, &row,&element,0.0,1.0,1.0e8); // large cost186 double element = 1.0; 187 build.addColumn(1, &row, &element, 0.0, 1.0, 1.0e8); // large cost 188 188 row++; 189 189 } 190 190 } 191 191 } 192 int nOK =0;193 int nNewRow =0;194 for (i =0;i<numberRows;i++) {192 int nOK = 0; 193 int nNewRow = 0; 194 for (i = 0; i < numberRows; i++) { 195 195 if (rowActivity[i]) 196 196 nOK++; 197 197 if (!rowActivity2[i]) 198 whichRow[nNewRow++] =i; // not satisfied198 whichRow[nNewRow++] = i; // not satisfied 199 199 else 200 modelPtr_>setRowStatus(i, ClpSimplex::basic); // make slack basic201 } 202 if (nOK <numberRows) {203 for (i =0;i<numberColumns;i++) {200 modelPtr_>setRowStatus(i, ClpSimplex::basic); // make slack basic 201 } 202 if (nOK < numberRows) { 203 for (i = 0; i < numberColumns; i++) { 204 204 if (!mark[i]) { 205 205 CoinBigIndex j; 206 int good =0;207 for (j =columnStart[i];208 j <columnStart[i]+columnLength[i];j++) {209 int iRow =row[j];206 int good = 0; 207 for (j = columnStart[i]; 208 j < columnStart[i] + columnLength[i]; j++) { 209 int iRow = row[j]; 210 210 if (!rowActivity[iRow]) { 211 rowActivity[iRow] 211 rowActivity[iRow]++; 212 212 good++; 213 213 } 214 214 } 215 215 if (good) { 216 nOK +=good;217 whichColumn[nNewCol++] =i;216 nOK += good; 217 whichColumn[nNewCol++] = i; 218 218 } 219 219 } 220 220 } 221 221 } 222 delete 223 delete 224 if (nOK <numberRows) {222 delete[] rowActivity; 223 delete[] rowActivity2; 224 if (nOK < numberRows) { 225 225 modelPtr_>setProblemStatus(1); 226 delete 227 delete 228 delete 226 delete[] whichRow; 227 delete[] whichColumn; 228 delete[] mark; 229 229 printf("infeasible by inspection\n"); 230 230 return; 231 231 } 232 bool allIn =false;233 if (nNewCol +nNoGood+numberRows>numberColumns) {232 bool allIn = false; 233 if (nNewCol + nNoGood + numberRows > numberColumns) { 234 234 // add in all 235 allIn =true;236 for (i =0;i<numberColumns;i++) {235 allIn = true; 236 for (i = 0; i < numberColumns; i++) { 237 237 if (!mark[i]) { 238 whichColumn[nNewCol++] =i;239 } 240 } 241 } 242 ClpSimplex * temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);238 whichColumn[nNewCol++] = i; 239 } 240 } 241 } 242 ClpSimplex *temp = new ClpSimplex(modelPtr_, nNewRow, whichRow, nNewCol, whichColumn); 243 243 if (equality) 244 244 temp>addColumns(build); 245 245 temp>setLogLevel(1); 246 246 printf("small has %d rows and %d columns (%d impossible to help) %s\n", 247 nNewRow,nNewCol,nNoGood,allIn ? "all in" : "");248 temp>setSpecialOptions(128 +512);247 nNewRow, nNewCol, nNoGood, allIn ? "all in" : ""); 248 temp>setSpecialOptions(128 + 512); 249 249 temp>setDualObjectiveLimit(1.0e50); 250 250 temp>dual(); 251 assert 252 double * 253 const double * 254 memset(solution, 0,numberColumns*sizeof(double));255 for (i =0;i<nNewCol;i++) {251 assert(!temp>status()); 252 double *solution = modelPtr_>primalColumnSolution(); 253 const double *solution2 = temp>primalColumnSolution(); 254 memset(solution, 0, numberColumns * sizeof(double)); 255 for (i = 0; i < nNewCol; i++) { 256 256 int iColumn = whichColumn[i]; 257 solution[iColumn] =solution2[i];258 modelPtr_>setStatus(iColumn, temp>getStatus(i));259 } 260 double * 261 const double * 262 double * 263 const double * 264 memset(dual, 0,numberRows*sizeof(double));265 for (i =0;i<nNewRow;i++) {266 int iRow =whichRow[i];267 modelPtr_>setRowStatus(iRow, temp>getRowStatus(i));268 rowSolution[iRow] =rowSolution2[i];269 dual[iRow] =dual2[i];257 solution[iColumn] = solution2[i]; 258 modelPtr_>setStatus(iColumn, temp>getStatus(i)); 259 } 260 double *rowSolution = modelPtr_>primalRowSolution(); 261 const double *rowSolution2 = temp>primalRowSolution(); 262 double *dual = modelPtr_>dualRowSolution(); 263 const double *dual2 = temp>dualRowSolution(); 264 memset(dual, 0, numberRows * sizeof(double)); 265 for (i = 0; i < nNewRow; i++) { 266 int iRow = whichRow[i]; 267 modelPtr_>setRowStatus(iRow, temp>getRowStatus(i)); 268 rowSolution[iRow] = rowSolution2[i]; 269 dual[iRow] = dual2[i]; 270 270 } 271 271 // See if optimal 272 double * 272 double *dj = modelPtr_>dualColumnSolution(); 273 273 // get reduced cost for large problem 274 274 // this assumes minimization 275 memcpy(dj, modelPtr_>objective(),numberColumns*sizeof(double));276 modelPtr_>transposeTimes(1.0, dual,dj);275 memcpy(dj, modelPtr_>objective(), numberColumns * sizeof(double)); 276 modelPtr_>transposeTimes(1.0, dual, dj); 277 277 modelPtr_>setObjectiveValue(temp>objectiveValue()); 278 278 modelPtr_>setProblemStatus(0); 279 int nBad =0;280 for (i =0;i<numberColumns;i++) {281 if (modelPtr_>getStatus(i) ==ClpSimplex::atLowerBound282 &&upper[i]>lower[i]&&dj[i]<1.0e5)279 int nBad = 0; 280 for (i = 0; i < numberColumns; i++) { 281 if (modelPtr_>getStatus(i) == ClpSimplex::atLowerBound 282 && upper[i] > lower[i] && dj[i] < 1.0e5) 283 283 nBad++; 284 284 } … … 286 286 //temp>writeMps("badb.mps"); 287 287 delete temp; 288 if (nBad &&!allIn) {289 assert (algorithm_==2);288 if (nBad && !allIn) { 289 assert(algorithm_ == 2); 290 290 //printf("%d bad\n",nBad); 291 291 timesBad_++; 292 292 // just non mark==2 293 int nAdded =0;294 for (i =0;i<numberColumns;i++) {293 int nAdded = 0; 294 for (i = 0; i < numberColumns; i++) { 295 295 if (!mark[i]) { 296 whichColumn[nNewCol++] =i;296 whichColumn[nNewCol++] = i; 297 297 nAdded++; 298 298 } 299 299 } 300 assert 300 assert(nAdded); 301 301 { 302 temp = new ClpSimplex(modelPtr_, nNewRow,whichRow,nNewCol,whichColumn);302 temp = new ClpSimplex(modelPtr_, nNewRow, whichRow, nNewCol, whichColumn); 303 303 if (equality) 304 304 temp>addColumns(build); 305 305 temp>setLogLevel(2); 306 temp>setSpecialOptions(128 +512);306 temp>setSpecialOptions(128 + 512); 307 307 temp>setDualObjectiveLimit(1.0e50); 308 308 temp>primal(1); 309 assert 310 double * 311 const double * 312 memset(solution, 0,numberColumns*sizeof(double));313 for (i =0;i<nNewCol;i++) {309 assert(!temp>status()); 310 double *solution = modelPtr_>primalColumnSolution(); 311 const double *solution2 = temp>primalColumnSolution(); 312 memset(solution, 0, numberColumns * sizeof(double)); 313 for (i = 0; i < nNewCol; i++) { 314 314 int iColumn = whichColumn[i]; 315 solution[iColumn] =solution2[i];316 modelPtr_>setStatus(iColumn, temp>getStatus(i));317 } 318 double * 319 const double * 320 double * 321 const double * 322 memset(dual, 0,numberRows*sizeof(double));323 for (i =0;i<nNewRow;i++) {324 int iRow =whichRow[i];325 modelPtr_>setRowStatus(iRow, temp>getRowStatus(i));326 rowSolution[iRow] =rowSolution2[i];327 dual[iRow] =dual2[i];315 solution[iColumn] = solution2[i]; 316 modelPtr_>setStatus(iColumn, temp>getStatus(i)); 317 } 318 double *rowSolution = modelPtr_>primalRowSolution(); 319 const double *rowSolution2 = temp>primalRowSolution(); 320 double *dual = modelPtr_>dualRowSolution(); 321 const double *dual2 = temp>dualRowSolution(); 322 memset(dual, 0, numberRows * sizeof(double)); 323 for (i = 0; i < nNewRow; i++) { 324 int iRow = whichRow[i]; 325 modelPtr_>setRowStatus(iRow, temp>getRowStatus(i)); 326 rowSolution[iRow] = rowSolution2[i]; 327 dual[iRow] = dual2[i]; 328 328 } 329 329 modelPtr_>setObjectiveValue(temp>objectiveValue()); 330 330 modelPtr_>setProblemStatus(0); 331 331 iterationsBad_ += temp>numberIterations(); 332 printf("clean %d\n", temp>numberIterations());332 printf("clean %d\n", temp>numberIterations()); 333 333 delete temp; 334 334 } 335 335 } 336 delete 337 delete 338 delete 336 delete[] mark; 337 delete[] whichRow; 338 delete[] whichColumn; 339 339 basis_ = getBasis(modelPtr_); 340 340 modelPtr_>setSpecialOptions(0); 341 341 count_++; 342 if ((count_ %100)==0&&algorithm_==2)343 printf("count %d, bad %d  iterations %d\n", count_,timesBad_,iterationsBad_);344 for (i =0;i<numberColumns;i++) {345 if (solution[i] >1.0e6modelPtr_>getStatus(i)==ClpSimplex::basic) {346 node_[i]=CoinMax(count_,node_[i]);347 348 } 349 } 350 if (modelPtr_>objectiveValue() >=modelPtr_>dualObjectiveLimit())342 if ((count_ % 100) == 0 && algorithm_ == 2) 343 printf("count %d, bad %d  iterations %d\n", count_, timesBad_, iterationsBad_); 344 for (i = 0; i < numberColumns; i++) { 345 if (solution[i] > 1.0e6  modelPtr_>getStatus(i) == ClpSimplex::basic) { 346 node_[i] = CoinMax(count_, node_[i]); 347 howMany_[i]++; 348 } 349 } 350 if (modelPtr_>objectiveValue() >= modelPtr_>dualObjectiveLimit()) 351 351 modelPtr_>setProblemStatus(1); 352 352 } … … 358 358 359 359 // 360 // Default Constructor 361 // 362 CbcSolver2::CbcSolver2 360 // Default Constructor 361 // 362 CbcSolver2::CbcSolver2() 363 363 : OsiClpSolverInterface() 364 364 { 365 node_ =NULL;366 howMany_ =NULL;367 count_ =0;365 node_ = NULL; 366 howMany_ = NULL; 367 count_ = 0; 368 368 model_ = NULL; 369 memory_ =300;370 algorithm_ =0;371 strategy_ =0;369 memory_ = 300; 370 algorithm_ = 0; 371 strategy_ = 0; 372 372 } 373 373 … … 375 375 // Clone 376 376 // 377 OsiSolverInterface * 377 OsiSolverInterface * 378 378 CbcSolver2::clone(bool CopyData) const 379 379 { … … 386 386 } 387 387 388 389 // 390 // Copy constructor 391 // 392 CbcSolver2::CbcSolver2 ( 393 const CbcSolver2 & rhs) 388 // 389 // Copy constructor 390 // 391 CbcSolver2::CbcSolver2( 392 const CbcSolver2 &rhs) 394 393 : OsiClpSolverInterface(rhs) 395 394 { 396 395 model_ = rhs.model_; 397 396 int numberColumns = modelPtr_>numberColumns(); 398 node_ =CoinCopyOfArray(rhs.node_,numberColumns);399 howMany_ =CoinCopyOfArray(rhs.howMany_,numberColumns);400 count_ =rhs.count_;401 memory_ =rhs.memory_;402 algorithm_ =rhs.algorithm_;403 strategy_ =rhs.strategy_;404 } 405 406 // 407 // Destructor 408 // 409 CbcSolver2::~CbcSolver2 410 { 411 delete 412 delete 413 } 414 415 // 416 // Assignment operator 397 node_ = CoinCopyOfArray(rhs.node_, numberColumns); 398 howMany_ = CoinCopyOfArray(rhs.howMany_, numberColumns); 399 count_ = rhs.count_; 400 memory_ = rhs.memory_; 401 algorithm_ = rhs.algorithm_; 402 strategy_ = rhs.strategy_; 403 } 404 405 // 406 // Destructor 407 // 408 CbcSolver2::~CbcSolver2() 409 { 410 delete[] node_; 411 delete[] howMany_; 412 } 413 414 // 415 // Assignment operator 417 416 // 418 417 CbcSolver2 & 419 CbcSolver2::operator=(const CbcSolver2 &rhs)420 { 421 if (this != &rhs) { 418 CbcSolver2::operator=(const CbcSolver2 &rhs) 419 { 420 if (this != &rhs) { 422 421 OsiClpSolverInterface::operator=(rhs); 423 delete 424 delete 422 delete[] node_; 423 delete[] howMany_; 425 424 model_ = rhs.model_; 426 425 int numberColumns = modelPtr_>numberColumns(); 427 node_ =CoinCopyOfArray(rhs.node_,numberColumns);428 howMany_ =CoinCopyOfArray(rhs.howMany_,numberColumns);429 count_ =rhs.count_;430 memory_ =rhs.memory_;431 algorithm_ =rhs.algorithm_;432 strategy_ =rhs.strategy_;426 node_ = CoinCopyOfArray(rhs.node_, numberColumns); 427 howMany_ = CoinCopyOfArray(rhs.howMany_, numberColumns); 428 count_ = rhs.count_; 429 memory_ = rhs.memory_; 430 algorithm_ = rhs.algorithm_; 431 strategy_ = rhs.strategy_; 433 432 } 434 433 return *this; … … 437 436 // Real initializer 438 437 // 439 void 440 CbcSolver2::initialize (CbcModel * model, const char * keep) 441 { 442 model_=model; 438 void CbcSolver2::initialize(CbcModel *model, const char *keep) 439 { 440 model_ = model; 443 441 int numberColumns = modelPtr_>numberColumns(); 444 442 if (numberColumns) { 445 443 node_ = new int[numberColumns]; 446 444 howMany_ = new int[numberColumns]; 447 for (int i =0;i<numberColumns;i++) {448 if (keep &&keep[i])449 node_[i]=COIN_INT_MAX;445 for (int i = 0; i < numberColumns; i++) { 446 if (keep && keep[i]) 447 node_[i] = COIN_INT_MAX; 450 448 else 451 node_[i]=0;452 howMany_[i] =0;449 node_[i] = 0; 450 howMany_[i] = 0; 453 451 } 454 452 } else { 455 node_ =NULL;456 howMany_ =NULL;453 node_ = NULL; 454 howMany_ = NULL; 457 455 } 458 456 }
Note: See TracChangeset
for help on using the changeset viewer.