Changeset 1286 for branches/sandbox/Cbc/src/CbcFathomDynamicProgramming.cpp
 Timestamp:
 Nov 9, 2009 6:33:07 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/sandbox/Cbc/src/CbcFathomDynamicProgramming.cpp
r1271 r1286 19 19 #include "CoinSort.hpp" 20 20 // Default Constructor 21 CbcFathomDynamicProgramming::CbcFathomDynamicProgramming() 22 :CbcFathom(),23 size_(0),24 type_(1),25 cost_(NULL),26 back_(NULL),27 lookup_(NULL),28 indices_(NULL),29 numberActive_(0),30 maximumSizeAllowed_(1000000),31 startBit_(NULL),32 numberBits_(NULL),33 rhs_(NULL),34 coefficients_(NULL),35 target_(0),36 numberNonOne_(0),37 bitPattern_(0),38 algorithm_(1)21 CbcFathomDynamicProgramming::CbcFathomDynamicProgramming() 22 : CbcFathom(), 23 size_(0), 24 type_(1), 25 cost_(NULL), 26 back_(NULL), 27 lookup_(NULL), 28 indices_(NULL), 29 numberActive_(0), 30 maximumSizeAllowed_(1000000), 31 startBit_(NULL), 32 numberBits_(NULL), 33 rhs_(NULL), 34 coefficients_(NULL), 35 target_(0), 36 numberNonOne_(0), 37 bitPattern_(0), 38 algorithm_(1) 39 39 { 40 40 … … 43 43 // Constructor from model 44 44 CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(CbcModel & model) 45 :CbcFathom(model),46 cost_(NULL),47 back_(NULL),48 lookup_(NULL),49 indices_(NULL),50 numberActive_(0),51 maximumSizeAllowed_(1000000),52 startBit_(NULL),53 numberBits_(NULL),54 rhs_(NULL),55 coefficients_(NULL),56 target_(0),57 numberNonOne_(0),58 bitPattern_(0),59 algorithm_(1)60 { 61 type_=checkPossible();62 } 63 64 // Destructor 45 : CbcFathom(model), 46 cost_(NULL), 47 back_(NULL), 48 lookup_(NULL), 49 indices_(NULL), 50 numberActive_(0), 51 maximumSizeAllowed_(1000000), 52 startBit_(NULL), 53 numberBits_(NULL), 54 rhs_(NULL), 55 coefficients_(NULL), 56 target_(0), 57 numberNonOne_(0), 58 bitPattern_(0), 59 algorithm_(1) 60 { 61 type_ = checkPossible(); 62 } 63 64 // Destructor 65 65 CbcFathomDynamicProgramming::~CbcFathomDynamicProgramming () 66 66 { 67 gutsOfDelete();67 gutsOfDelete(); 68 68 } 69 69 // Does deleteions 70 void 70 void 71 71 CbcFathomDynamicProgramming::gutsOfDelete() 72 72 { 73 delete [] cost_;74 delete [] back_;75 delete [] lookup_;76 delete [] indices_;77 delete [] startBit_;78 delete [] numberBits_;79 delete [] rhs_;80 delete [] coefficients_;81 cost_ = NULL;82 back_ = NULL;83 lookup_ = NULL;84 indices_ =NULL;85 startBit_ = NULL;86 numberBits_ = NULL;87 rhs_ = NULL;88 coefficients_ = NULL;73 delete [] cost_; 74 delete [] back_; 75 delete [] lookup_; 76 delete [] indices_; 77 delete [] startBit_; 78 delete [] numberBits_; 79 delete [] rhs_; 80 delete [] coefficients_; 81 cost_ = NULL; 82 back_ = NULL; 83 lookup_ = NULL; 84 indices_ = NULL; 85 startBit_ = NULL; 86 numberBits_ = NULL; 87 rhs_ = NULL; 88 coefficients_ = NULL; 89 89 } 90 90 // Clone … … 92 92 CbcFathomDynamicProgramming::clone() const 93 93 { 94 return new CbcFathomDynamicProgramming(*this);95 } 96 97 // Copy constructor 94 return new CbcFathomDynamicProgramming(*this); 95 } 96 97 // Copy constructor 98 98 CbcFathomDynamicProgramming::CbcFathomDynamicProgramming(const CbcFathomDynamicProgramming & rhs) 99 :100 CbcFathom(rhs),101 size_(rhs.size_),102 type_(rhs.type_),103 cost_(NULL),104 back_(NULL),105 lookup_(NULL),106 indices_(NULL),107 numberActive_(rhs.numberActive_),108 maximumSizeAllowed_(rhs.maximumSizeAllowed_),109 startBit_(NULL),110 numberBits_(NULL),111 rhs_(NULL),112 coefficients_(NULL),113 target_(rhs.target_),114 numberNonOne_(rhs.numberNonOne_),115 bitPattern_(rhs.bitPattern_),116 algorithm_(rhs.algorithm_)117 { 118 if (size_) {119 cost_=CoinCopyOfArray(rhs.cost_,size_);120 back_=CoinCopyOfArray(rhs.back_,size_);121 int numberRows=model_>getNumRows();122 lookup_=CoinCopyOfArray(rhs.lookup_,numberRows);123 startBit_=CoinCopyOfArray(rhs.startBit_,numberActive_);124 indices_=CoinCopyOfArray(rhs.indices_,numberActive_);125 numberBits_=CoinCopyOfArray(rhs.numberBits_,numberActive_);126 rhs_=CoinCopyOfArray(rhs.rhs_,numberActive_);127 coefficients_=CoinCopyOfArray(rhs.coefficients_,numberActive_);128 }99 : 100 CbcFathom(rhs), 101 size_(rhs.size_), 102 type_(rhs.type_), 103 cost_(NULL), 104 back_(NULL), 105 lookup_(NULL), 106 indices_(NULL), 107 numberActive_(rhs.numberActive_), 108 maximumSizeAllowed_(rhs.maximumSizeAllowed_), 109 startBit_(NULL), 110 numberBits_(NULL), 111 rhs_(NULL), 112 coefficients_(NULL), 113 target_(rhs.target_), 114 numberNonOne_(rhs.numberNonOne_), 115 bitPattern_(rhs.bitPattern_), 116 algorithm_(rhs.algorithm_) 117 { 118 if (size_) { 119 cost_ = CoinCopyOfArray(rhs.cost_, size_); 120 back_ = CoinCopyOfArray(rhs.back_, size_); 121 int numberRows = model_>getNumRows(); 122 lookup_ = CoinCopyOfArray(rhs.lookup_, numberRows); 123 startBit_ = CoinCopyOfArray(rhs.startBit_, numberActive_); 124 indices_ = CoinCopyOfArray(rhs.indices_, numberActive_); 125 numberBits_ = CoinCopyOfArray(rhs.numberBits_, numberActive_); 126 rhs_ = CoinCopyOfArray(rhs.rhs_, numberActive_); 127 coefficients_ = CoinCopyOfArray(rhs.coefficients_, numberActive_); 128 } 129 129 } 130 130 // Returns type 131 int 131 int 132 132 CbcFathomDynamicProgramming::checkPossible(int allowableSize) 133 133 { 134 algorithm_=1; 135 assert(model_>solver()); 136 OsiSolverInterface * solver = model_>solver(); 137 const CoinPackedMatrix * matrix = solver>getMatrixByCol(); 138 139 int numberIntegers = model_>numberIntegers(); 140 int numberColumns = solver>getNumCols(); 141 size_=0; 142 if (numberIntegers!=numberColumns) 143 return 1; // can't do dynamic programming 144 145 const double * lower = solver>getColLower(); 146 const double * upper = solver>getColUpper(); 147 const double * rowUpper = solver>getRowUpper(); 148 149 int numberRows = model_>getNumRows(); 150 int i; 151 152 // First check columns to see if possible 153 double * rhs = new double [numberRows]; 154 CoinCopyN(rowUpper,numberRows,rhs); 155 156 // Column copy 157 const double * element = matrix>getElements(); 158 const int * row = matrix>getIndices(); 159 const CoinBigIndex * columnStart = matrix>getVectorStarts(); 160 const int * columnLength = matrix>getVectorLengths(); 161 bool bad=false; 162 /* It is just possible that we could say okay as 163 variables may get fixed but seems unlikely */ 164 for (i=0;i<numberColumns;i++) { 165 int j; 166 double lowerValue = lower[i]; 167 assert (lowerValue==floor(lowerValue)); 168 for (j=columnStart[i]; 169 j<columnStart[i]+columnLength[i];j++) { 170 int iRow=row[j]; 171 double value = element[j]; 172 if (upper[i]>lowerValue&&(value<=0.0value!=floor(value))) 173 bad=true; 174 if (lowerValue) 175 rhs[iRow] = lowerValue*value; 176 } 177 } 178 // check possible (at present do not allow covering) 179 int numberActive=0; 180 bool infeasible=false; 181 bool saveBad=bad; 182 for (i=0;i<numberRows;i++) { 183 if (rhs[i]<0) 184 infeasible=true; 185 else if (rhs[i]>1.0e5fabs(rhs[i]floor(rhs[i]+0.5))>1.0e7) 186 bad=true; 187 else if (rhs[i]>0.0) 188 numberActive++; 189 } 190 if (badinfeasible) { 191 delete [] rhs; 192 if (!saveBad&&infeasible) 193 return 2; 194 else 195 return 1; 196 } 197 // check size of array needed 198 double size=1.0; 199 double check = COIN_INT_MAX; 200 for (i=0;i<numberRows;i++) { 201 int n= static_cast<int> (floor(rhs[i]+0.5)); 202 if (n) { 203 n++; // allow for 0,1... n 204 if (numberActive!=1) { 205 // power of 2 206 int iBit=0; 207 int k=n; 208 k &= ~1; 209 while (k) { 210 iBit++; 211 k &= ~(1<<iBit); 212 } 213 // See if exact power 214 if (n!=(1<<iBit)) { 215 // round up to next power of 2 216 n= 1<<(iBit+1); 217 } 218 size *= n; 219 if (size>=check) 220 break; 221 } else { 222 size = n; // just one constraint 223 } 224 } 225 } 226 // set size needed 227 if (size>=check) 228 size_=COIN_INT_MAX; 229 else 230 size_=static_cast<int> (size); 231 232 int n01=0; 233 int nbadcoeff=0; 234 // See if we can tighten bounds 235 for (i=0;i<numberColumns;i++) { 236 int j; 237 double lowerValue = lower[i]; 238 double gap = upper[i]lowerValue; 239 for (j=columnStart[i]; 240 j<columnStart[i]+columnLength[i];j++) { 241 int iRow=row[j]; 242 double value = element[j]; 243 if (value!=1.0) 244 nbadcoeff++; 245 if (gap*value>rhs[iRow]+1.0e8) 246 gap = rhs[iRow]/value; 247 } 248 gap=lowerValue+floor(gap+1.0e7); 249 if (gap<upper[i]) 250 solver>setColUpper(i,gap); 251 if (gap<=1.0) 252 n01++; 253 } 254 if (allowableSize&&size_<=allowableSize) { 255 if (n01==numberColumns&&!nbadcoeff) 256 algorithm_ = 0; // easiest 257 else 258 algorithm_ = 1; 259 } 260 if (allowableSize&&size_<=allowableSize) { 261 numberActive_=numberActive; 262 indices_ = new int [numberActive_]; 263 cost_ = new double [size_]; 264 CoinFillN(cost_,size_,COIN_DBL_MAX); 265 // but do nothing is okay 266 cost_[0]=0.0; 267 back_ = new int[size_]; 268 CoinFillN(back_,size_,1); 269 startBit_=new int[numberActive_]; 270 numberBits_=new int[numberActive_]; 271 lookup_ = new int [numberRows]; 272 rhs_ = new int [numberActive_]; 273 numberActive=0; 274 int kBit=0; 275 for (i=0;i<numberRows;i++) { 276 int n= static_cast<int> (floor(rhs[i]+0.5)); 277 if (n) { 278 lookup_[i]=numberActive; 279 rhs_[numberActive]=n; 280 startBit_[numberActive]=kBit; 281 n++; // allow for 0,1... n 282 int iBit=0; 283 // power of 2 284 int k=n; 285 k &= ~1; 286 while (k) { 287 iBit++; 288 k &= ~(1<<iBit); 289 } 290 // See if exact power 291 if (n!=(1<<iBit)) { 292 // round up to next power of 2 293 iBit++; 294 } 295 if (numberActive!=1) { 296 n= 1<<iBit; 297 size *= n; 298 if (size>=check) 299 break; 300 } else { 301 size = n; // just one constraint 302 } 303 numberBits_[numberActive++]=iBit; 304 kBit += iBit; 305 } else { 306 lookup_[i]=1; 307 } 308 } 309 const double * rowLower = solver>getRowLower(); 310 if (algorithm_==0) { 311 // rhs 1 and coefficients 1 312 // Get first possible solution for printing 313 target_=1; 314 int needed=0; 315 int numberActive=0; 316 for (i=0;i<numberRows;i++) { 317 int newRow = lookup_[i]; 318 if (newRow>=0) { 319 if (rowLower[i]==rowUpper[i]) { 320 needed += 1<<numberActive; 321 numberActive++; 322 } 323 } 324 } 325 for (i=0;i<size_;i++) { 326 if ((i&needed)==needed) { 327 break; 328 } 329 } 330 target_=i; 331 } else { 332 coefficients_ = new int[numberActive_]; 333 // If not too many general rhs then we can be more efficient 334 numberNonOne_=0; 335 for (i=0;i<numberActive_;i++) { 336 if (rhs_[i]!=1) 337 numberNonOne_++; 338 } 339 if (numberNonOne_*2<numberActive_) { 340 // put rhs >1 every second 341 int * permute = new int[numberActive_]; 342 int * temp = new int[numberActive_]; 343 // try different ways 344 int k=0; 345 for (i=0;i<numberRows;i++) { 346 int newRow = lookup_[i]; 347 if (newRow>=0&&rhs_[newRow]>1) { 348 permute[newRow]=k; 349 k +=2; 350 } 351 } 352 // adjust so k points to last 353 k = 2; 354 // and now rest 355 int k1=1; 356 for (i=0;i<numberRows;i++) { 357 int newRow = lookup_[i]; 358 if (newRow>=0&&rhs_[newRow]==1) { 359 permute[newRow]=k1; 360 k1++; 361 if (k1<=k) 362 k1++; 363 } 364 } 365 for (i=0;i<numberActive_;i++) { 366 int put = permute[i]; 367 temp[put]=rhs_[i]; 368 } 369 memcpy(rhs_,temp,numberActive_*sizeof(int)); 370 for (i=0;i<numberActive_;i++) { 371 int put = permute[i]; 372 temp[put]=numberBits_[i]; 373 } 374 memcpy(numberBits_,temp,numberActive_*sizeof(int)); 375 k=0; 376 for (i=0;i<numberActive_;i++) { 377 startBit_[i]=k; 378 k+= numberBits_[i]; 379 } 380 for (i=0;i<numberRows;i++) { 381 int newRow = lookup_[i]; 382 if (newRow>=0) 383 lookup_[i]=permute[newRow]; 384 } 385 delete [] permute; 386 delete [] temp; 387 // mark new method 388 algorithm_=2; 389 } 390 // Get first possible solution for printing 391 target_=1; 392 int needed=0; 393 int * lower2 = new int[numberActive_]; 394 for (i=0;i<numberRows;i++) { 395 int newRow = lookup_[i]; 396 if (newRow>=0) { 397 int gap=static_cast<int> (rowUpper[i]CoinMax(0.0,rowLower[i])); 398 lower2[newRow]=rhs_[newRow]gap; 399 int numberBits = numberBits_[newRow]; 400 int startBit = startBit_[newRow]; 401 if (numberBits==1&&!gap) { 402 needed = 1<<startBit; 403 } 404 } 405 } 406 for (i=0;i<size_;i++) { 407 if ((i&needed)==needed) { 408 // this one may do 409 bool good = true; 410 for (int kk=0;kk<numberActive_;kk++) { 411 int numberBits = numberBits_[kk]; 412 int startBit = startBit_[kk]; 413 int size = 1<<numberBits; 414 int start = 1<<startBit; 415 int mask = start*(size1); 416 int level=(i&mask)>>startBit; 417 if (level<lower2[kk]) { 418 good=false; 419 break; 420 } 421 } 422 if (good) { 423 break; 424 } 425 } 426 } 427 delete [] lower2; 428 target_=i; 429 } 430 } 431 delete [] rhs; 432 if (allowableSize&&size_>allowableSize) { 433 printf("Too large  need %d entries x 8 bytes\n",size_); 434 return 1; // too big 435 } else { 436 return algorithm_; 437 } 438 } 439 440 // Resets stuff if model changes 441 void 442 CbcFathomDynamicProgramming::resetModel(CbcModel * model) 443 { 444 model_=model; 445 type_=checkPossible(); 446 } 447 int 448 CbcFathomDynamicProgramming::fathom(double * & betterSolution) 449 { 450 int returnCode=0; 451 int type=checkPossible(maximumSizeAllowed_); 452 assert (type!=1); 453 if (type==2) { 454 // infeasible (so complete search done) 455 return 1; 456 } 457 if (algorithm_>=0) { 134 algorithm_ = 1; 135 assert(model_>solver()); 458 136 OsiSolverInterface * solver = model_>solver(); 137 const CoinPackedMatrix * matrix = solver>getMatrixByCol(); 138 139 int numberIntegers = model_>numberIntegers(); 140 int numberColumns = solver>getNumCols(); 141 size_ = 0; 142 if (numberIntegers != numberColumns) 143 return 1; // can't do dynamic programming 144 459 145 const double * lower = solver>getColLower(); 460 146 const double * upper = solver>getColUpper(); 461 const double * objective = solver>getObjCoefficients(); 462 double direction = solver>getObjSense(); 463 const CoinPackedMatrix * matrix = solver>getMatrixByCol(); 147 const double * rowUpper = solver>getRowUpper(); 148 149 int numberRows = model_>getNumRows(); 150 int i; 151 152 // First check columns to see if possible 153 double * rhs = new double [numberRows]; 154 CoinCopyN(rowUpper, numberRows, rhs); 155 464 156 // Column copy 465 157 const double * element = matrix>getElements(); … … 467 159 const CoinBigIndex * columnStart = matrix>getVectorStarts(); 468 160 const int * columnLength = matrix>getVectorLengths(); 469 const double * rowLower = solver>getRowLower(); 470 const double * rowUpper = solver>getRowUpper(); 471 int numberRows = model_>getNumRows(); 472 473 int numberColumns = solver>getNumCols(); 474 double offset; 475 solver>getDblParam(OsiObjOffset,offset); 476 double fixedObj=offset; 477 int i; 478 // may be possible 479 double bestAtTarget = COIN_DBL_MAX; 480 for (i=0;i<numberColumns;i++) { 481 if (size_>10000000&&(i%100)==0) 482 printf("column %d\n",i); 483 double lowerValue = lower[i]; 484 assert (lowerValue==floor(lowerValue)); 485 double cost = direction * objective[i]; 486 fixedObj += lowerValue*cost; 487 int gap = static_cast<int> (upper[i]lowerValue); 488 CoinBigIndex start = columnStart[i]; 489 tryColumn(columnLength[i],row+start,element+start,cost,gap); 490 if (cost_[target_]<bestAtTarget) { 491 if (model_>messageHandler()>logLevel()>1) 492 printf("At column %d new best objective of %g\n",i,cost_[target_]); 493 bestAtTarget = cost_[target_]; 494 } 495 } 496 returnCode=1; 497 int needed=0; 498 double bestValue=COIN_DBL_MAX; 499 int iBest=1; 500 if (algorithm_==0) { 501 int numberActive=0; 502 for (i=0;i<numberRows;i++) { 503 int newRow = lookup_[i]; 504 if (newRow>=0) { 505 if (rowLower[i]==rowUpper[i]) { 506 needed += 1<<numberActive; 161 bool bad = false; 162 /* It is just possible that we could say okay as 163 variables may get fixed but seems unlikely */ 164 for (i = 0; i < numberColumns; i++) { 165 int j; 166 double lowerValue = lower[i]; 167 assert (lowerValue == floor(lowerValue)); 168 for (j = columnStart[i]; 169 j < columnStart[i] + columnLength[i]; j++) { 170 int iRow = row[j]; 171 double value = element[j]; 172 if (upper[i] > lowerValue && (value <= 0.0  value != floor(value))) 173 bad = true; 174 if (lowerValue) 175 rhs[iRow] = lowerValue * value; 176 } 177 } 178 // check possible (at present do not allow covering) 179 int numberActive = 0; 180 bool infeasible = false; 181 bool saveBad = bad; 182 for (i = 0; i < numberRows; i++) { 183 if (rhs[i] < 0) 184 infeasible = true; 185 else if (rhs[i] > 1.0e5  fabs(rhs[i]  floor(rhs[i] + 0.5)) > 1.0e7) 186 bad = true; 187 else if (rhs[i] > 0.0) 507 188 numberActive++; 508 } 509 } 510 } 511 for (i=0;i<size_;i++) { 512 if ((i&needed)==needed) { 513 // this one will do 514 if (cost_[i]<bestValue) { 515 bestValue=cost_[i]; 516 iBest=i; 517 } 518 } 519 } 189 } 190 if (bad  infeasible) { 191 delete [] rhs; 192 if (!saveBad && infeasible) 193 return 2; 194 else 195 return 1; 196 } 197 // check size of array needed 198 double size = 1.0; 199 double check = COIN_INT_MAX; 200 for (i = 0; i < numberRows; i++) { 201 int n = static_cast<int> (floor(rhs[i] + 0.5)); 202 if (n) { 203 n++; // allow for 0,1... n 204 if (numberActive != 1) { 205 // power of 2 206 int iBit = 0; 207 int k = n; 208 k &= ~1; 209 while (k) { 210 iBit++; 211 k &= ~(1 << iBit); 212 } 213 // See if exact power 214 if (n != (1 << iBit)) { 215 // round up to next power of 2 216 n = 1 << (iBit + 1); 217 } 218 size *= n; 219 if (size >= check) 220 break; 221 } else { 222 size = n; // just one constraint 223 } 224 } 225 } 226 // set size needed 227 if (size >= check) 228 size_ = COIN_INT_MAX; 229 else 230 size_ = static_cast<int> (size); 231 232 int n01 = 0; 233 int nbadcoeff = 0; 234 // See if we can tighten bounds 235 for (i = 0; i < numberColumns; i++) { 236 int j; 237 double lowerValue = lower[i]; 238 double gap = upper[i]  lowerValue; 239 for (j = columnStart[i]; 240 j < columnStart[i] + columnLength[i]; j++) { 241 int iRow = row[j]; 242 double value = element[j]; 243 if (value != 1.0) 244 nbadcoeff++; 245 if (gap*value > rhs[iRow] + 1.0e8) 246 gap = rhs[iRow] / value; 247 } 248 gap = lowerValue + floor(gap + 1.0e7); 249 if (gap < upper[i]) 250 solver>setColUpper(i, gap); 251 if (gap <= 1.0) 252 n01++; 253 } 254 if (allowableSize && size_ <= allowableSize) { 255 if (n01 == numberColumns && !nbadcoeff) 256 algorithm_ = 0; // easiest 257 else 258 algorithm_ = 1; 259 } 260 if (allowableSize && size_ <= allowableSize) { 261 numberActive_ = numberActive; 262 indices_ = new int [numberActive_]; 263 cost_ = new double [size_]; 264 CoinFillN(cost_, size_, COIN_DBL_MAX); 265 // but do nothing is okay 266 cost_[0] = 0.0; 267 back_ = new int[size_]; 268 CoinFillN(back_, size_, 1); 269 startBit_ = new int[numberActive_]; 270 numberBits_ = new int[numberActive_]; 271 lookup_ = new int [numberRows]; 272 rhs_ = new int [numberActive_]; 273 numberActive = 0; 274 int kBit = 0; 275 for (i = 0; i < numberRows; i++) { 276 int n = static_cast<int> (floor(rhs[i] + 0.5)); 277 if (n) { 278 lookup_[i] = numberActive; 279 rhs_[numberActive] = n; 280 startBit_[numberActive] = kBit; 281 n++; // allow for 0,1... n 282 int iBit = 0; 283 // power of 2 284 int k = n; 285 k &= ~1; 286 while (k) { 287 iBit++; 288 k &= ~(1 << iBit); 289 } 290 // See if exact power 291 if (n != (1 << iBit)) { 292 // round up to next power of 2 293 iBit++; 294 } 295 if (numberActive != 1) { 296 n = 1 << iBit; 297 size *= n; 298 if (size >= check) 299 break; 300 } else { 301 size = n; // just one constraint 302 } 303 numberBits_[numberActive++] = iBit; 304 kBit += iBit; 305 } else { 306 lookup_[i] = 1; 307 } 308 } 309 const double * rowLower = solver>getRowLower(); 310 if (algorithm_ == 0) { 311 // rhs 1 and coefficients 1 312 // Get first possible solution for printing 313 target_ = 1; 314 int needed = 0; 315 int numberActive = 0; 316 for (i = 0; i < numberRows; i++) { 317 int newRow = lookup_[i]; 318 if (newRow >= 0) { 319 if (rowLower[i] == rowUpper[i]) { 320 needed += 1 << numberActive; 321 numberActive++; 322 } 323 } 324 } 325 for (i = 0; i < size_; i++) { 326 if ((i&needed) == needed) { 327 break; 328 } 329 } 330 target_ = i; 331 } else { 332 coefficients_ = new int[numberActive_]; 333 // If not too many general rhs then we can be more efficient 334 numberNonOne_ = 0; 335 for (i = 0; i < numberActive_; i++) { 336 if (rhs_[i] != 1) 337 numberNonOne_++; 338 } 339 if (numberNonOne_*2 < numberActive_) { 340 // put rhs >1 every second 341 int * permute = new int[numberActive_]; 342 int * temp = new int[numberActive_]; 343 // try different ways 344 int k = 0; 345 for (i = 0; i < numberRows; i++) { 346 int newRow = lookup_[i]; 347 if (newRow >= 0 && rhs_[newRow] > 1) { 348 permute[newRow] = k; 349 k += 2; 350 } 351 } 352 // adjust so k points to last 353 k = 2; 354 // and now rest 355 int k1 = 1; 356 for (i = 0; i < numberRows; i++) { 357 int newRow = lookup_[i]; 358 if (newRow >= 0 && rhs_[newRow] == 1) { 359 permute[newRow] = k1; 360 k1++; 361 if (k1 <= k) 362 k1++; 363 } 364 } 365 for (i = 0; i < numberActive_; i++) { 366 int put = permute[i]; 367 temp[put] = rhs_[i]; 368 } 369 memcpy(rhs_, temp, numberActive_*sizeof(int)); 370 for (i = 0; i < numberActive_; i++) { 371 int put = permute[i]; 372 temp[put] = numberBits_[i]; 373 } 374 memcpy(numberBits_, temp, numberActive_*sizeof(int)); 375 k = 0; 376 for (i = 0; i < numberActive_; i++) { 377 startBit_[i] = k; 378 k += numberBits_[i]; 379 } 380 for (i = 0; i < numberRows; i++) { 381 int newRow = lookup_[i]; 382 if (newRow >= 0) 383 lookup_[i] = permute[newRow]; 384 } 385 delete [] permute; 386 delete [] temp; 387 // mark new method 388 algorithm_ = 2; 389 } 390 // Get first possible solution for printing 391 target_ = 1; 392 int needed = 0; 393 int * lower2 = new int[numberActive_]; 394 for (i = 0; i < numberRows; i++) { 395 int newRow = lookup_[i]; 396 if (newRow >= 0) { 397 int gap = static_cast<int> (rowUpper[i]  CoinMax(0.0, rowLower[i])); 398 lower2[newRow] = rhs_[newRow]  gap; 399 int numberBits = numberBits_[newRow]; 400 int startBit = startBit_[newRow]; 401 if (numberBits == 1 && !gap) { 402 needed = 1 << startBit; 403 } 404 } 405 } 406 for (i = 0; i < size_; i++) { 407 if ((i&needed) == needed) { 408 // this one may do 409 bool good = true; 410 for (int kk = 0; kk < numberActive_; kk++) { 411 int numberBits = numberBits_[kk]; 412 int startBit = startBit_[kk]; 413 int size = 1 << numberBits; 414 int start = 1 << startBit; 415 int mask = start * (size  1); 416 int level = (i & mask) >> startBit; 417 if (level < lower2[kk]) { 418 good = false; 419 break; 420 } 421 } 422 if (good) { 423 break; 424 } 425 } 426 } 427 delete [] lower2; 428 target_ = i; 429 } 430 } 431 delete [] rhs; 432 if (allowableSize && size_ > allowableSize) { 433 printf("Too large  need %d entries x 8 bytes\n", size_); 434 return 1; // too big 520 435 } else { 521 int * lower = new int[numberActive_]; 522 for (i=0;i<numberRows;i++) { 523 int newRow = lookup_[i]; 524 if (newRow>=0) { 525 int gap=static_cast<int> (rowUpper[i]CoinMax(0.0,rowLower[i])); 526 lower[newRow]=rhs_[newRow]gap; 527 int numberBits = numberBits_[newRow]; 528 int startBit = startBit_[newRow]; 529 if (numberBits==1&&!gap) { 530 needed = 1<<startBit; 531 } 532 } 533 } 534 for (i=0;i<size_;i++) { 535 if ((i&needed)==needed) { 536 // this one may do 537 bool good = true; 538 for (int kk=0;kk<numberActive_;kk++) { 539 int numberBits = numberBits_[kk]; 540 int startBit = startBit_[kk]; 541 int size = 1<<numberBits; 542 int start = 1<<startBit; 543 int mask = start*(size1); 544 int level=(i&mask)>>startBit; 545 if (level<lower[kk]) { 546 good=false; 547 break; 548 } 549 } 550 if (good&&cost_[i]<bestValue) { 551 bestValue=cost_[i]; 552 iBest=i; 553 } 554 } 555 } 556 delete [] lower; 557 } 558 if (bestValue<COIN_DBL_MAX) { 559 bestValue += fixedObj; 560 if (model_>messageHandler()>logLevel()>1) 561 printf("Can get solution of %g\n",bestValue); 562 if (bestValue<model_>getMinimizationObjValue()) { 563 // set up solution 564 betterSolution = new double[numberColumns]; 565 memcpy(betterSolution,lower,numberColumns*sizeof(double)); 566 while (iBest>0) { 567 int n=decodeBitPattern(iBestback_[iBest],indices_,numberRows); 568 // Search for cheapest 569 double bestCost=COIN_DBL_MAX; 570 int iColumn=1; 571 for (i=0;i<numberColumns;i++) { 572 if (n==columnLength[i]) { 573 bool good=true; 574 for (int j=columnStart[i]; 575 j<columnStart[i]+columnLength[i];j++) { 576 int iRow=row[j]; 577 double value = element[j]; 578 int iValue = static_cast<int> (value); 579 if(iValue!=indices_[iRow]) { 580 good=false; 581 break; 582 } 583 } 584 if (good && objective[i]<bestCost&&betterSolution[i]<upper[i]) { 585 bestCost=objective[i]; 586 iColumn=i; 587 } 588 } 589 } 590 assert (iColumn>=0); 591 betterSolution[iColumn]++; 592 assert (betterSolution[iColumn]<=upper[iColumn]); 593 iBest = back_[iBest]; 594 } 595 } 596 // paranoid check 597 double * rowActivity = new double [numberRows]; 598 memset(rowActivity,0,numberRows*sizeof(double)); 599 for (i=0;i<numberColumns;i++) { 600 int j; 601 double value = betterSolution[i]; 602 if (value) { 603 for (j=columnStart[i]; 604 j<columnStart[i]+columnLength[i];j++) { 605 int iRow=row[j]; 606 rowActivity[iRow] += value*element[j]; 607 } 608 } 609 } 610 // check was feasible 611 bool feasible=true; 612 for (i=0;i<numberRows;i++) { 613 if(rowActivity[i]<rowLower[i]) { 614 if (rowActivity[i]<rowLower[i]1.0e8) 615 feasible = false; 616 } else if(rowActivity[i]>rowUpper[i]) { 617 if (rowActivity[i]>rowUpper[i]+1.0e8) 618 feasible = false; 619 } 620 } 621 if (feasible) { 622 if (model_>messageHandler()>logLevel()>0) 623 printf("** good solution of %g by dynamic programming\n",bestValue); 624 } 625 delete [] rowActivity; 626 } 627 gutsOfDelete(); 628 } 629 return returnCode; 436 return algorithm_; 437 } 438 } 439 440 // Resets stuff if model changes 441 void 442 CbcFathomDynamicProgramming::resetModel(CbcModel * model) 443 { 444 model_ = model; 445 type_ = checkPossible(); 446 } 447 int 448 CbcFathomDynamicProgramming::fathom(double * & betterSolution) 449 { 450 int returnCode = 0; 451 int type = checkPossible(maximumSizeAllowed_); 452 assert (type != 1); 453 if (type == 2) { 454 // infeasible (so complete search done) 455 return 1; 456 } 457 if (algorithm_ >= 0) { 458 OsiSolverInterface * solver = model_>solver(); 459 const double * lower = solver>getColLower(); 460 const double * upper = solver>getColUpper(); 461 const double * objective = solver>getObjCoefficients(); 462 double direction = solver>getObjSense(); 463 const CoinPackedMatrix * matrix = solver>getMatrixByCol(); 464 // Column copy 465 const double * element = matrix>getElements(); 466 const int * row = matrix>getIndices(); 467 const CoinBigIndex * columnStart = matrix>getVectorStarts(); 468 const int * columnLength = matrix>getVectorLengths(); 469 const double * rowLower = solver>getRowLower(); 470 const double * rowUpper = solver>getRowUpper(); 471 int numberRows = model_>getNumRows(); 472 473 int numberColumns = solver>getNumCols(); 474 double offset; 475 solver>getDblParam(OsiObjOffset, offset); 476 double fixedObj = offset; 477 int i; 478 // may be possible 479 double bestAtTarget = COIN_DBL_MAX; 480 for (i = 0; i < numberColumns; i++) { 481 if (size_ > 10000000 && (i % 100) == 0) 482 printf("column %d\n", i); 483 double lowerValue = lower[i]; 484 assert (lowerValue == floor(lowerValue)); 485 double cost = direction * objective[i]; 486 fixedObj += lowerValue * cost; 487 int gap = static_cast<int> (upper[i]  lowerValue); 488 CoinBigIndex start = columnStart[i]; 489 tryColumn(columnLength[i], row + start, element + start, cost, gap); 490 if (cost_[target_] < bestAtTarget) { 491 if (model_>messageHandler()>logLevel() > 1) 492 printf("At column %d new best objective of %g\n", i, cost_[target_]); 493 bestAtTarget = cost_[target_]; 494 } 495 } 496 returnCode = 1; 497 int needed = 0; 498 double bestValue = COIN_DBL_MAX; 499 int iBest = 1; 500 if (algorithm_ == 0) { 501 int numberActive = 0; 502 for (i = 0; i < numberRows; i++) { 503 int newRow = lookup_[i]; 504 if (newRow >= 0) { 505 if (rowLower[i] == rowUpper[i]) { 506 needed += 1 << numberActive; 507 numberActive++; 508 } 509 } 510 } 511 for (i = 0; i < size_; i++) { 512 if ((i&needed) == needed) { 513 // this one will do 514 if (cost_[i] < bestValue) { 515 bestValue = cost_[i]; 516 iBest = i; 517 } 518 } 519 } 520 } else { 521 int * lower = new int[numberActive_]; 522 for (i = 0; i < numberRows; i++) { 523 int newRow = lookup_[i]; 524 if (newRow >= 0) { 525 int gap = static_cast<int> (rowUpper[i]  CoinMax(0.0, rowLower[i])); 526 lower[newRow] = rhs_[newRow]  gap; 527 int numberBits = numberBits_[newRow]; 528 int startBit = startBit_[newRow]; 529 if (numberBits == 1 && !gap) { 530 needed = 1 << startBit; 531 } 532 } 533 } 534 for (i = 0; i < size_; i++) { 535 if ((i&needed) == needed) { 536 // this one may do 537 bool good = true; 538 for (int kk = 0; kk < numberActive_; kk++) { 539 int numberBits = numberBits_[kk]; 540 int startBit = startBit_[kk]; 541 int size = 1 << numberBits; 542 int start = 1 << startBit; 543 int mask = start * (size  1); 544 int level = (i & mask) >> startBit; 545 if (level < lower[kk]) { 546 good = false; 547 break; 548 } 549 } 550 if (good && cost_[i] < bestValue) { 551 bestValue = cost_[i]; 552 iBest = i; 553 } 554 } 555 } 556 delete [] lower; 557 } 558 if (bestValue < COIN_DBL_MAX) { 559 bestValue += fixedObj; 560 if (model_>messageHandler()>logLevel() > 1) 561 printf("Can get solution of %g\n", bestValue); 562 if (bestValue < model_>getMinimizationObjValue()) { 563 // set up solution 564 betterSolution = new double[numberColumns]; 565 memcpy(betterSolution, lower, numberColumns*sizeof(double)); 566 while (iBest > 0) { 567 int n = decodeBitPattern(iBest  back_[iBest], indices_, numberRows); 568 // Search for cheapest 569 double bestCost = COIN_DBL_MAX; 570 int iColumn = 1; 571 for (i = 0; i < numberColumns; i++) { 572 if (n == columnLength[i]) { 573 bool good = true; 574 for (int j = columnStart[i]; 575 j < columnStart[i] + columnLength[i]; j++) { 576 int iRow = row[j]; 577 double value = element[j]; 578 int iValue = static_cast<int> (value); 579 if (iValue != indices_[iRow]) { 580 good = false; 581 break; 582 } 583 } 584 if (good && objective[i] < bestCost && betterSolution[i] < upper[i]) { 585 bestCost = objective[i]; 586 iColumn = i; 587 } 588 } 589 } 590 assert (iColumn >= 0); 591 betterSolution[iColumn]++; 592 assert (betterSolution[iColumn] <= upper[iColumn]); 593 iBest = back_[iBest]; 594 } 595 } 596 // paranoid check 597 double * rowActivity = new double [numberRows]; 598 memset(rowActivity, 0, numberRows*sizeof(double)); 599 for (i = 0; i < numberColumns; i++) { 600 int j; 601 double value = betterSolution[i]; 602 if (value) { 603 for (j = columnStart[i]; 604 j < columnStart[i] + columnLength[i]; j++) { 605 int iRow = row[j]; 606 rowActivity[iRow] += value * element[j]; 607 } 608 } 609 } 610 // check was feasible 611 bool feasible = true; 612 for (i = 0; i < numberRows; i++) { 613 if (rowActivity[i] < rowLower[i]) { 614 if (rowActivity[i] < rowLower[i]  1.0e8) 615 feasible = false; 616 } else if (rowActivity[i] > rowUpper[i]) { 617 if (rowActivity[i] > rowUpper[i] + 1.0e8) 618 feasible = false; 619 } 620 } 621 if (feasible) { 622 if (model_>messageHandler()>logLevel() > 0) 623 printf("** good solution of %g by dynamic programming\n", bestValue); 624 } 625 delete [] rowActivity; 626 } 627 gutsOfDelete(); 628 } 629 return returnCode; 630 630 } 631 631 /* Tries a column 632 632 returns true if was used in making any changes. 633 633 */ 634 bool 634 bool 635 635 CbcFathomDynamicProgramming::tryColumn(int numberElements, const int * rows, 636 637 638 { 639 bool touched=false;640 int n=0;641 if (algorithm_==0) {642 for (int j=0;j<numberElements;j++) {643 int iRow=rows[j];644 double value = coefficients[j];645 int newRow = lookup_[iRow];646 if (newRow<0value>rhs_[newRow]) {647 n=0;648 break; //can't use649 } else {650 indices_[n++]=newRow;651 }652 }653 if (n&&upper) {654 touched=addOneColumn0(n,indices_,cost);655 }656 } else {657 for (int j=0;j<numberElements;j++) {658 int iRow=rows[j];659 double value = coefficients[j];660 int iValue = static_cast<int> (value);661 int newRow = lookup_[iRow];662 if (newRow<0iValue>rhs_[newRow]) {663 n=0;664 break; //can't use665 } else {666 coefficients_[n]=iValue;667 indices_[n++]=newRow;668 if (upper*iValue>rhs_[newRow]) {669 upper = rhs_[newRow]/iValue;670 }671 }672 }673 if (n) {674 if (algorithm_==1) {675 for (int k=1;k<=upper;k++) {676 bool t = addOneColumn1(n,indices_,coefficients_,cost);677 if (t)678 touched=true;679 }680 } else {681 CoinSort_2(indices_,indices_+n,coefficients_);682 for (int k=1;k<=upper;k++) {683 bool t = addOneColumn1A(n,indices_,coefficients_,cost);684 if (t)685 touched=true;686 }687 }688 }689 }690 return touched;636 const double * coefficients, double cost, 637 int upper) 638 { 639 bool touched = false; 640 int n = 0; 641 if (algorithm_ == 0) { 642 for (int j = 0; j < numberElements; j++) { 643 int iRow = rows[j]; 644 double value = coefficients[j]; 645 int newRow = lookup_[iRow]; 646 if (newRow < 0  value > rhs_[newRow]) { 647 n = 0; 648 break; //can't use 649 } else { 650 indices_[n++] = newRow; 651 } 652 } 653 if (n && upper) { 654 touched = addOneColumn0(n, indices_, cost); 655 } 656 } else { 657 for (int j = 0; j < numberElements; j++) { 658 int iRow = rows[j]; 659 double value = coefficients[j]; 660 int iValue = static_cast<int> (value); 661 int newRow = lookup_[iRow]; 662 if (newRow < 0  iValue > rhs_[newRow]) { 663 n = 0; 664 break; //can't use 665 } else { 666 coefficients_[n] = iValue; 667 indices_[n++] = newRow; 668 if (upper*iValue > rhs_[newRow]) { 669 upper = rhs_[newRow] / iValue; 670 } 671 } 672 } 673 if (n) { 674 if (algorithm_ == 1) { 675 for (int k = 1; k <= upper; k++) { 676 bool t = addOneColumn1(n, indices_, coefficients_, cost); 677 if (t) 678 touched = true; 679 } 680 } else { 681 CoinSort_2(indices_, indices_ + n, coefficients_); 682 for (int k = 1; k <= upper; k++) { 683 bool t = addOneColumn1A(n, indices_, coefficients_, cost); 684 if (t) 685 touched = true; 686 } 687 } 688 } 689 } 690 return touched; 691 691 } 692 692 /* Adds one column if type 0, 693 693 returns true if was used in making any changes 694 694 */ 695 bool 695 bool 696 696 CbcFathomDynamicProgramming::addOneColumn0(int numberElements, const int * rows, 697 698 { 699 // build up mask700 int mask=0;701 int i;702 for (i=0;i<numberElements;i++) {703 int iRow=rows[i];704 mask = 1<<iRow;705 }706 bitPattern_=mask;707 i=size_1mask;708 bool touched = false;709 while (i>=0) {710 int kMask = i&mask;711 if (kMask==0) {712 double thisCost = cost_[i];713 if (thisCost!=COIN_DBL_MAX) {714 // possible715 double newCost=thisCost+cost;716 int next = i + mask;717 if (cost_[next]>newCost) {718 cost_[next]=newCost;719 back_[next]=i;720 touched=true;721 }722 }723 i;724 } else {725 // we can skip some726 int k=(i&~mask);697 double cost) 698 { 699 // build up mask 700 int mask = 0; 701 int i; 702 for (i = 0; i < numberElements; i++) { 703 int iRow = rows[i]; 704 mask = 1 << iRow; 705 } 706 bitPattern_ = mask; 707 i = size_  1  mask; 708 bool touched = false; 709 while (i >= 0) { 710 int kMask = i & mask; 711 if (kMask == 0) { 712 double thisCost = cost_[i]; 713 if (thisCost != COIN_DBL_MAX) { 714 // possible 715 double newCost = thisCost + cost; 716 int next = i + mask; 717 if (cost_[next] > newCost) { 718 cost_[next] = newCost; 719 back_[next] = i; 720 touched = true; 721 } 722 } 723 i; 724 } else { 725 // we can skip some 726 int k = (i&~mask); 727 727 #ifdef CBC_DEBUG 728 for (int j=i1;j>k;j) {729 int jMask = j&mask;730 assert (jMask!=0);731 }728 for (int j = i  1; j > k; j) { 729 int jMask = j & mask; 730 assert (jMask != 0); 731 } 732 732 #endif 733 i=k;734 }735 }736 return touched;733 i = k; 734 } 735 } 736 return touched; 737 737 } 738 738 /* Adds one attempt of one column of type 1, … … 740 740 At present the user has to call it once for each possible value 741 741 */ 742 bool 742 bool 743 743 CbcFathomDynamicProgramming::addOneColumn1(int numberElements, const int * rows, 744 745 { 746 /* build up masks.747 a) mask for 1 rhs748 b) mask for addition749 c) mask so adding will overflow750 d) individual masks751 */752 int mask1=0;753 int maskAdd=0;754 int mask2=0;755 int i;756 int n2=0;757 int mask[40];758 int nextbit[40];759 int adjust[40];760 assert (numberElements<=40);761 for (i=0;i<numberElements;i++) {762 int iRow=rows[i];763 int numberBits = numberBits_[iRow];764 int startBit = startBit_[iRow];765 if (numberBits==1) {766 mask1 = 1<<startBit;767 maskAdd = 1<<startBit;768 mask2 = 1<<startBit;769 } else {770 int value=coefficients[i];771 int size = 1<<numberBits;772 int start = 1<<startBit;773 assert (value<size);774 maskAdd = start*value;775 int gap = sizerhs_[iRow]1;776 assert (gap>=0);777 int hi2=rhs_[iRow]value;778 if (hi2<size1)779 hi2++;780 adjust[n2] = start*hi2;781 mask2 += start*gap;782 nextbit[n2]=startBit+numberBits;783 mask[n2++] = start*(size1);784 }785 }786 bitPattern_=maskAdd;787 i=size_1maskAdd;788 bool touched = false;789 while (i>=0) {790 int kMask = i&mask1;791 if (kMask==0) {792 bool good=true;793 for (int kk=n21;kk>=0;kk) {794 int iMask = mask[kk];795 int jMask = iMask&mask2;796 int kkMask = iMask&i;797 kkMask += jMask;798 if (kkMask>iMask) {799 // we can skip some800 int k=(i&~iMask);801 k = adjust[kk];744 const int * coefficients, double cost) 745 { 746 /* build up masks. 747 a) mask for 1 rhs 748 b) mask for addition 749 c) mask so adding will overflow 750 d) individual masks 751 */ 752 int mask1 = 0; 753 int maskAdd = 0; 754 int mask2 = 0; 755 int i; 756 int n2 = 0; 757 int mask[40]; 758 int nextbit[40]; 759 int adjust[40]; 760 assert (numberElements <= 40); 761 for (i = 0; i < numberElements; i++) { 762 int iRow = rows[i]; 763 int numberBits = numberBits_[iRow]; 764 int startBit = startBit_[iRow]; 765 if (numberBits == 1) { 766 mask1 = 1 << startBit; 767 maskAdd = 1 << startBit; 768 mask2 = 1 << startBit; 769 } else { 770 int value = coefficients[i]; 771 int size = 1 << numberBits; 772 int start = 1 << startBit; 773 assert (value < size); 774 maskAdd = start * value; 775 int gap = size  rhs_[iRow]  1; 776 assert (gap >= 0); 777 int hi2 = rhs_[iRow]  value; 778 if (hi2 < size  1) 779 hi2++; 780 adjust[n2] = start * hi2; 781 mask2 += start * gap; 782 nextbit[n2] = startBit + numberBits; 783 mask[n2++] = start * (size  1); 784 } 785 } 786 bitPattern_ = maskAdd; 787 i = size_  1  maskAdd; 788 bool touched = false; 789 while (i >= 0) { 790 int kMask = i & mask1; 791 if (kMask == 0) { 792 bool good = true; 793 for (int kk = n2  1; kk >= 0; kk) { 794 int iMask = mask[kk]; 795 int jMask = iMask & mask2; 796 int kkMask = iMask & i; 797 kkMask += jMask; 798 if (kkMask > iMask) { 799 // we can skip some 800 int k = (i&~iMask); 801 k = adjust[kk]; 802 802 #ifdef CBC_DEBUG 803 for (int j=i1;j>k;j) {804 int jMask = j&mask1;805 if (jMask==0) {806 bool good=true;807 for (int kk=n21;kk>=0;kk) {808 int iMask = mask[kk];809 int jMask = iMask&mask2;810 int kkMask = iMask&i;811 kkMask += jMask;812 if (kkMask>iMask) {813 good=false;814 break;815 }816 }817 assert (!good);818 }819 }803 for (int j = i  1; j > k; j) { 804 int jMask = j & mask1; 805 if (jMask == 0) { 806 bool good = true; 807 for (int kk = n2  1; kk >= 0; kk) { 808 int iMask = mask[kk]; 809 int jMask = iMask & mask2; 810 int kkMask = iMask & i; 811 kkMask += jMask; 812 if (kkMask > iMask) { 813 good = false; 814 break; 815 } 816 } 817 assert (!good); 818 } 819 } 820 820 #endif 821 i=k;822 good=false;823 break;824 }825 }826 if (good) {827 double thisCost = cost_[i];828 if (thisCost!=COIN_DBL_MAX) {829 // possible830 double newCost=thisCost+cost;831 int next = i + maskAdd;832 if (cost_[next]>newCost) {833 cost_[next]=newCost;834 back_[next]=i;835 touched=true;836 }837 }838 }839 i;840 } else {841 // we can skip some842 // we can skip some843 int k=(i&~mask1);821 i = k; 822 good = false; 823 break; 824 } 825 } 826 if (good) { 827 double thisCost = cost_[i]; 828 if (thisCost != COIN_DBL_MAX) { 829 // possible 830 double newCost = thisCost + cost; 831 int next = i + maskAdd; 832 if (cost_[next] > newCost) { 833 cost_[next] = newCost; 834 back_[next] = i; 835 touched = true; 836 } 837 } 838 } 839 i; 840 } else { 841 // we can skip some 842 // we can skip some 843 int k = (i&~mask1); 844 844 #ifdef CBC_DEBUG 845 for (int j=i1;j>k;j) {846 int jMask = j&mask1;847 assert (jMask!=0);848 }845 for (int j = i  1; j > k; j) { 846 int jMask = j & mask1; 847 assert (jMask != 0); 848 } 849 849 #endif 850 i=k;851 }852 }853 return touched;850 i = k; 851 } 852 } 853 return touched; 854 854 } 855 855 /* Adds one attempt of one column of type 1, … … 858 858 This version is when there are enough 1 rhs to do faster 859 859 */ 860 bool 860 bool 861 861 CbcFathomDynamicProgramming::addOneColumn1A(int numberElements, const int * rows, 862 const int * coefficients, double cost) 863 { 864 /* build up masks. 865 a) mask for 1 rhs 866 b) mask for addition 867 c) mask so adding will overflow 868 d) mask for non 1 rhs 869 */ 870 int maskA=0; 871 int maskAdd=0; 872 int maskC=0; 873 int maskD=0; 874 int i; 875 for (i=0;i<numberElements;i++) { 876 int iRow=rows[i]; 877 int numberBits = numberBits_[iRow]; 878 int startBit = startBit_[iRow]; 879 if (numberBits==1) { 880 maskA = 1<<startBit; 881 maskAdd = 1<<startBit; 862 const int * coefficients, double cost) 863 { 864 /* build up masks. 865 a) mask for 1 rhs 866 b) mask for addition 867 c) mask so adding will overflow 868 d) mask for non 1 rhs 869 */ 870 int maskA = 0; 871 int maskAdd = 0; 872 int maskC = 0; 873 int maskD = 0; 874 int i; 875 for (i = 0; i < numberElements; i++) { 876 int iRow = rows[i]; 877 int numberBits = numberBits_[iRow]; 878 int startBit = startBit_[iRow]; 879 if (numberBits == 1) { 880 maskA = 1 << startBit; 881 maskAdd = 1 << startBit; 882 } else { 883 int value = coefficients[i]; 884 int size = 1 << numberBits; 885 int start = 1 << startBit; 886 assert (value < size); 887 maskAdd = start * value; 888 int gap = size  rhs_[iRow] + value  1; 889 assert (gap > 0 && gap <= size  1); 890 maskC = start * gap; 891 maskD = start * (size  1); 892 } 893 } 894 bitPattern_ = maskAdd; 895 int maskDiff = maskD  maskC; 896 i = size_  1  maskAdd; 897 bool touched = false; 898 if (!maskD) { 899 // Just ones 900 while (i >= 0) { 901 int kMask = i & maskA; 902 if (kMask == 0) { 903 double thisCost = cost_[i]; 904 if (thisCost != COIN_DBL_MAX) { 905 // possible 906 double newCost = thisCost + cost; 907 int next = i + maskAdd; 908 if (cost_[next] > newCost) { 909 cost_[next] = newCost; 910 back_[next] = i; 911 touched = true; 912 } 913 } 914 i; 915 } else { 916 // we can skip some 917 int k = (i&~maskA); 918 i = k; 919 } 920 } 882 921 } else { 883 int value=coefficients[i]; 884 int size = 1<<numberBits; 885 int start = 1<<startBit; 886 assert (value<size); 887 maskAdd = start*value; 888 int gap = sizerhs_[iRow]+value1; 889 assert (gap>0&&gap<=size1); 890 maskC = start*gap; 891 maskD = start*(size1); 892 } 893 } 894 bitPattern_=maskAdd; 895 int maskDiff = maskDmaskC; 896 i=size_1maskAdd; 897 bool touched = false; 898 if (!maskD) { 899 // Just ones 900 while (i>=0) { 901 int kMask = i&maskA; 902 if (kMask==0) { 903 double thisCost = cost_[i]; 904 if (thisCost!=COIN_DBL_MAX) { 905 // possible 906 double newCost=thisCost+cost; 907 int next = i + maskAdd; 908 if (cost_[next]>newCost) { 909 cost_[next]=newCost; 910 back_[next]=i; 911 touched=true; 912 } 913 } 914 i; 915 } else { 916 // we can skip some 917 int k=(i&~maskA); 918 i=k; 919 } 920 } 921 } else { 922 // More general 923 while (i>=0) { 924 int kMask = i&maskA; 925 if (kMask==0) { 926 int added = i & maskD; // just bits belonging to non 1 rhs 927 added += maskC; // will overflow mask if bad 928 added &= (~maskD); 929 if (added == 0) { 930 double thisCost = cost_[i]; 931 if (thisCost!=COIN_DBL_MAX) { 932 // possible 933 double newCost=thisCost+cost; 934 int next = i + maskAdd; 935 if (cost_[next]>newCost) { 936 cost_[next]=newCost; 937 back_[next]=i; 938 touched=true; 939 } 940 } 941 i; 942 } else { 943 // we can skip some 944 int k = i & ~ maskD; // clear all 945 // Put back enough  but only below where we are 946 int kk=(numberNonOne_<<1)2; 947 assert (rhs_[kk]>1); 948 int iMask=0; 949 for(;kk>=0;kk=2) { 950 iMask = 1<<startBit_[kk+1]; 951 if ((added&iMask)!=0) { 952 iMask; 953 break; 954 } 955 } 956 assert (kk>=0); 957 iMask &= maskDiff; 958 k = iMask; 959 assert (k<i); 960 i=k; 961 } 962 } else { 963 // we can skip some 964 int k=(i&~maskA); 965 i=k; 966 } 967 } 968 } 969 return touched; 922 // More general 923 while (i >= 0) { 924 int kMask = i & maskA; 925 if (kMask == 0) { 926 int added = i & maskD; // just bits belonging to non 1 rhs 927 added += maskC; // will overflow mask if bad 928 added &= (~maskD); 929 if (added == 0) { 930 double thisCost = cost_[i]; 931 if (thisCost != COIN_DBL_MAX) { 932 // possible 933 double newCost = thisCost + cost; 934 int next = i + maskAdd; 935 if (cost_[next] > newCost) { 936 cost_[next] = newCost; 937 back_[next] = i; 938 touched = true; 939 } 940 } 941 i; 942 } else { 943 // we can skip some 944 int k = i & ~ maskD; // clear all 945 // Put back enough  but only below where we are 946 int kk = (numberNonOne_ << 1)  2; 947 assert (rhs_[kk] > 1); 948 int iMask = 0; 949 for (; kk >= 0; kk = 2) { 950 iMask = 1 << startBit_[kk+1]; 951 if ((added&iMask) != 0) { 952 iMask; 953 break; 954 } 955 } 956 assert (kk >= 0); 957 iMask &= maskDiff; 958 k = iMask; 959 assert (k < i); 960 i = k; 961 } 962 } else { 963 // we can skip some 964 int k = (i&~maskA); 965 i = k; 966 } 967 } 968 } 969 return touched; 970 970 } 971 971 // update model 972 972 void CbcFathomDynamicProgramming::setModel(CbcModel * model) 973 973 { 974 model_ = model;975 type_=checkPossible();974 model_ = model; 975 type_ = checkPossible(); 976 976 } 977 977 // Gets bit pattern from original column 978 978 int CbcFathomDynamicProgramming::bitPattern(int numberElements, const int * rows, 979 980 { 981 int i;982 int mask=0;983 switch (algorithm_) {984 // just ones985 case 0:986 for (i=0;i<numberElements;i++) {987 int iRow=rows[i];988 iRow = lookup_[iRow];989 if (iRow>=0)990 mask = 1<<iRow;991 }992 break;993 //994 case 1:995 case 2:996 for (i=0;i<numberElements;i++) {997 int iRow=rows[i];998 iRow = lookup_[iRow];999 if (iRow>=0) {1000 int startBit = startBit_[iRow];1001 int value=coefficients[i];1002 int start = 1<<startBit;1003 mask = start*value;1004 }1005 }1006 break;1007 }1008 return mask;979 const int * coefficients) 980 { 981 int i; 982 int mask = 0; 983 switch (algorithm_) { 984 // just ones 985 case 0: 986 for (i = 0; i < numberElements; i++) { 987 int iRow = rows[i]; 988 iRow = lookup_[iRow]; 989 if (iRow >= 0) 990 mask = 1 << iRow; 991 } 992 break; 993 // 994 case 1: 995 case 2: 996 for (i = 0; i < numberElements; i++) { 997 int iRow = rows[i]; 998 iRow = lookup_[iRow]; 999 if (iRow >= 0) { 1000 int startBit = startBit_[iRow]; 1001 int value = coefficients[i]; 1002 int start = 1 << startBit; 1003 mask = start * value; 1004 } 1005 } 1006 break; 1007 } 1008 return mask; 1009 1009 } 1010 1010 // Fills in original column (dense) from bit pattern 1011 1011 int CbcFathomDynamicProgramming::decodeBitPattern(int bitPattern, 1012 1013 1014 { 1015 int i;1016 int n=0;1017 switch (algorithm_) {1018 // just ones1019 case 0:1020 for (i=0;i<numberRows;i++) {1021 values[i]=0;1022 int iRow = lookup_[i];1023 if (iRow>=0) {1024 if ((bitPattern&(1<<iRow))!=0) {1025 values[i]=1;1026 n++;1027 }1028 }1029 }1030 break;1031 //1032 case 1:1033 case 2:1034 for (i=0;i<numberRows;i++) {1035 values[i]=0;1036 int iRow = lookup_[i];1037 if (iRow>=0) {1038 int startBit = startBit_[iRow];1039 int numberBits = numberBits_[iRow];1040 int iValue = bitPattern>>startBit;1041 iValue &= ((1<<numberBits)1);1042 if (iValue) {1043 values[i]=iValue;1044 n++;1045 }1046 }1047 }1048 break;1049 }1050 return n;1051 } 1052 1053 1012 int * values, 1013 int numberRows) 1014 { 1015 int i; 1016 int n = 0; 1017 switch (algorithm_) { 1018 // just ones 1019 case 0: 1020 for (i = 0; i < numberRows; i++) { 1021 values[i] = 0; 1022 int iRow = lookup_[i]; 1023 if (iRow >= 0) { 1024 if ((bitPattern&(1 << iRow)) != 0) { 1025 values[i] = 1; 1026 n++; 1027 } 1028 } 1029 } 1030 break; 1031 // 1032 case 1: 1033 case 2: 1034 for (i = 0; i < numberRows; i++) { 1035 values[i] = 0; 1036 int iRow = lookup_[i]; 1037 if (iRow >= 0) { 1038 int startBit = startBit_[iRow]; 1039 int numberBits = numberBits_[iRow]; 1040 int iValue = bitPattern >> startBit; 1041 iValue &= ((1 << numberBits)  1); 1042 if (iValue) { 1043 values[i] = iValue; 1044 n++; 1045 } 1046 } 1047 } 1048 break; 1049 } 1050 return n; 1051 } 1052 1053
Note: See TracChangeset
for help on using the changeset viewer.