Changeset 425 for branches/devel/Cbc/examples/CbcBranchLink.cpp
 Timestamp:
 Sep 15, 2006 4:57:25 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel/Cbc/examples/CbcBranchLink.cpp
r424 r425 15 15 #include "CbcBranchLink.hpp" 16 16 #include "CoinError.hpp" 17 #include "CoinPackedMatrix.hpp" 17 18 18 19 // Default Constructor … … 22 23 numberMembers_(0), 23 24 numberLinks_(0), 24 which_(NULL) 25 which_(NULL), 26 sosType_(1) 25 27 { 26 28 } … … 32 34 numberMembers_(numberMembers), 33 35 numberLinks_(numberLinks), 34 which_(NULL) 36 which_(NULL), 37 sosType_(1) 35 38 { 36 39 id_=identifier; … … 61 64 // Useful constructor (which are indices) 62 65 CbcLink::CbcLink (CbcModel * model, int numberMembers, 63 int numberLinks, const int * which , const double * weights, int identifier)66 int numberLinks, int sosType, const int * which , const double * weights, int identifier) 64 67 : CbcObject(model), 65 68 numberMembers_(numberMembers), 66 69 numberLinks_(numberLinks), 67 which_(NULL) 70 which_(NULL), 71 sosType_(sosType) 68 72 { 69 73 id_=identifier; … … 98 102 numberMembers_ = rhs.numberMembers_; 99 103 numberLinks_ = rhs.numberLinks_; 104 sosType_ = rhs.sosType_; 100 105 if (numberMembers_) { 101 106 weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_); … … 124 129 numberMembers_ = rhs.numberMembers_; 125 130 numberLinks_ = rhs.numberLinks_; 131 sosType_ = rhs.sosType_; 126 132 if (numberMembers_) { 127 133 weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_); … … 180 186 #endif 181 187 } 182 weight += weights_[j]*CoinMin(value,upper[iColumn]); 188 value = CoinMin(value,upper[iColumn]); 189 weight += weights_[j]*value; 183 190 if (firstNonZero<0) 184 191 firstNonZero=j; … … 188 195 base += numberLinks_; 189 196 } 197 double valueInfeasibility; 190 198 preferredWay=1; 191 if (lastNonZerofirstNonZero>= 1) {199 if (lastNonZerofirstNonZero>=sosType_) { 192 200 // find where to branch 193 201 assert (sum>0.0); 194 202 weight /= sum; 195 double value = lastNonZerofirstNonZero+1; 196 value *= 0.5/((double) numberMembers_); 197 return value; 203 valueInfeasibility = lastNonZerofirstNonZero+1; 204 valueInfeasibility *= 0.5/((double) numberMembers_); 205 //#define DISTANCE 206 #ifdef DISTANCE 207 assert (sosType_==1); // code up 208 /* may still be satisfied. 209 For LOS type 2 we might wish to move coding around 210 and keep initial info in model_ for speed 211 */ 212 int iWhere; 213 bool possible=false; 214 for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) { 215 if (fabs(weightweights_[iWhere])<1.0e8) { 216 possible=true; 217 break; 218 } 219 } 220 if (possible) { 221 // One could move some of this (+ arrays) into model_ 222 const CoinPackedMatrix * matrix = solver>getMatrixByCol(); 223 const double * element = matrix>getMutableElements(); 224 const int * row = matrix>getIndices(); 225 const CoinBigIndex * columnStart = matrix>getVectorStarts(); 226 const int * columnLength = matrix>getVectorLengths(); 227 const double * rowSolution = solver>getRowActivity(); 228 const double * rowLower = solver>getRowLower(); 229 const double * rowUpper = solver>getRowUpper(); 230 int numberRows = matrix>getNumRows(); 231 double * array = new double [numberRows]; 232 CoinZeroN(array,numberRows); 233 int * which = new int [numberRows]; 234 int n=0; 235 int base=numberLinks_*firstNonZero; 236 for (j=firstNonZero;j<=lastNonZero;j++) { 237 for (int k=0;k<numberLinks_;k++) { 238 int iColumn = which_[base+k]; 239 double value = CoinMax(0.0,solution[iColumn]); 240 if (value>integerTolerance&&upper[iColumn]) { 241 value = CoinMin(value,upper[iColumn]); 242 for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { 243 int iRow = row[j]; 244 double a = array[iRow]; 245 if (a) { 246 a += value*element[j]; 247 if (!a) 248 a = 1.0e100; 249 } else { 250 which[n++]=iRow; 251 a=value*element[j]; 252 assert (a); 253 } 254 array[iRow]=a; 255 } 256 } 257 } 258 base += numberLinks_; 259 } 260 base=numberLinks_*iWhere; 261 for (int k=0;k<numberLinks_;k++) { 262 int iColumn = which_[base+k]; 263 const double value = 1.0; 264 for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { 265 int iRow = row[j]; 266 double a = array[iRow]; 267 if (a) { 268 a = value*element[j]; 269 if (!a) 270 a = 1.0e100; 271 } else { 272 which[n++]=iRow; 273 a=value*element[j]; 274 assert (a); 275 } 276 array[iRow]=a; 277 } 278 } 279 for (j=0;j<n;j++) { 280 int iRow = which[j]; 281 // moving to point will increase row solution by this 282 double distance = array[iRow]; 283 if (distance>1.0e8) { 284 if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e8) { 285 possible=false; 286 break; 287 } 288 } else if (distance<1.0e8) { 289 if (distance+rowSolution[iRow]<rowLower[iRow]1.0e8) { 290 possible=false; 291 break; 292 } 293 } 294 } 295 for (j=0;j<n;j++) 296 array[which[j]]=0.0; 297 delete [] array; 298 delete [] which; 299 if (possible) { 300 valueInfeasibility=0.0; 301 printf("possible %d %d %d\n",firstNonZero,lastNonZero,iWhere); 302 } 303 } 304 #endif 198 305 } else { 199 return 0.0; // satisfied 200 } 306 valueInfeasibility = 0.0; // satisfied 307 } 308 return valueInfeasibility; 201 309 } 202 310 … … 231 339 base += numberLinks_; 232 340 } 233 assert (lastNonZerofirstNonZero==0) ; 341 #ifdef DISTANCE 342 if (lastNonZerofirstNonZero>sosType_1) { 343 /* may still be satisfied. 344 For LOS type 2 we might wish to move coding around 345 and keep initial info in model_ for speed 346 */ 347 int iWhere; 348 bool possible=false; 349 for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) { 350 if (fabs(weightweights_[iWhere])<1.0e8) { 351 possible=true; 352 break; 353 } 354 } 355 if (possible) { 356 // One could move some of this (+ arrays) into model_ 357 const CoinPackedMatrix * matrix = solver>getMatrixByCol(); 358 const double * element = matrix>getMutableElements(); 359 const int * row = matrix>getIndices(); 360 const CoinBigIndex * columnStart = matrix>getVectorStarts(); 361 const int * columnLength = matrix>getVectorLengths(); 362 const double * rowSolution = solver>getRowActivity(); 363 const double * rowLower = solver>getRowLower(); 364 const double * rowUpper = solver>getRowUpper(); 365 int numberRows = matrix>getNumRows(); 366 double * array = new double [numberRows]; 367 CoinZeroN(array,numberRows); 368 int * which = new int [numberRows]; 369 int n=0; 370 int base=numberLinks_*firstNonZero; 371 for (j=firstNonZero;j<=lastNonZero;j++) { 372 for (int k=0;k<numberLinks_;k++) { 373 int iColumn = which_[base+k]; 374 double value = CoinMax(0.0,solution[iColumn]); 375 if (value>integerTolerance&&upper[iColumn]) { 376 value = CoinMin(value,upper[iColumn]); 377 for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { 378 int iRow = row[j]; 379 double a = array[iRow]; 380 if (a) { 381 a += value*element[j]; 382 if (!a) 383 a = 1.0e100; 384 } else { 385 which[n++]=iRow; 386 a=value*element[j]; 387 assert (a); 388 } 389 array[iRow]=a; 390 } 391 } 392 } 393 base += numberLinks_; 394 } 395 base=numberLinks_*iWhere; 396 for (int k=0;k<numberLinks_;k++) { 397 int iColumn = which_[base+k]; 398 const double value = 1.0; 399 for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { 400 int iRow = row[j]; 401 double a = array[iRow]; 402 if (a) { 403 a = value*element[j]; 404 if (!a) 405 a = 1.0e100; 406 } else { 407 which[n++]=iRow; 408 a=value*element[j]; 409 assert (a); 410 } 411 array[iRow]=a; 412 } 413 } 414 for (j=0;j<n;j++) { 415 int iRow = which[j]; 416 // moving to point will increase row solution by this 417 double distance = array[iRow]; 418 if (distance>1.0e8) { 419 if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e8) { 420 possible=false; 421 break; 422 } 423 } else if (distance<1.0e8) { 424 if (distance+rowSolution[iRow]<rowLower[iRow]1.0e8) { 425 possible=false; 426 break; 427 } 428 } 429 } 430 for (j=0;j<n;j++) 431 array[which[j]]=0.0; 432 delete [] array; 433 delete [] which; 434 if (possible) { 435 printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere); 436 firstNonZero=iWhere; 437 lastNonZero=iWhere; 438 } 439 } 440 } 441 #else 442 assert (lastNonZerofirstNonZero<sosType_) ; 443 #endif 234 444 base=0; 235 445 for (j=0;j<firstNonZero;j++) { … … 288 498 base += numberLinks_; 289 499 } 290 assert (lastNonZerofirstNonZero>= 1) ;500 assert (lastNonZerofirstNonZero>=sosType_) ; 291 501 // find where to branch 292 502 assert (sum>0.0); … … 297 507 if (weight<weights_[iWhere+1]) 298 508 break; 299 separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]); 509 if (sosType_==1) { 510 // SOS 1 511 separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]); 512 } else { 513 // SOS 2 514 if (iWhere==firstNonFixed) 515 iWhere++;; 516 if (iWhere==lastNonFixed1) 517 iWhere = lastNonFixed2; 518 separator = weights_[iWhere+1]; 519 } 300 520 // create object 301 521 CbcBranchingObject * branch;
Note: See TracChangeset
for help on using the changeset viewer.