[1574] | 1 | // $Id: CbcBranchLink.cpp 1902 2013-04-10 16:58:16Z stefan $ |
---|
[121] | 2 | // Copyright (C) 2005, International Business Machines |
---|
| 3 | // Corporation and others. All Rights Reserved. |
---|
[1574] | 4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
| 5 | |
---|
[121] | 6 | #include <cassert> |
---|
| 7 | #include <cmath> |
---|
| 8 | #include <cfloat> |
---|
| 9 | //#define CBC_DEBUG |
---|
| 10 | |
---|
[1902] | 11 | #include "CoinPragma.hpp" |
---|
[121] | 12 | #include "OsiSolverInterface.hpp" |
---|
| 13 | #include "CbcModel.hpp" |
---|
| 14 | #include "CbcMessage.hpp" |
---|
| 15 | #include "CbcBranchLink.hpp" |
---|
| 16 | #include "CoinError.hpp" |
---|
[640] | 17 | #include "CoinPackedMatrix.hpp" |
---|
[121] | 18 | |
---|
| 19 | // Default Constructor |
---|
| 20 | CbcLink::CbcLink () |
---|
| 21 | : CbcObject(), |
---|
| 22 | weights_(NULL), |
---|
| 23 | numberMembers_(0), |
---|
| 24 | numberLinks_(0), |
---|
[640] | 25 | which_(NULL), |
---|
| 26 | sosType_(1) |
---|
[121] | 27 | { |
---|
| 28 | } |
---|
| 29 | |
---|
| 30 | // Useful constructor (which are indices) |
---|
| 31 | CbcLink::CbcLink (CbcModel * model, int numberMembers, |
---|
| 32 | int numberLinks, int first , const double * weights, int identifier) |
---|
| 33 | : CbcObject(model), |
---|
| 34 | numberMembers_(numberMembers), |
---|
| 35 | numberLinks_(numberLinks), |
---|
[640] | 36 | which_(NULL), |
---|
| 37 | sosType_(1) |
---|
[121] | 38 | { |
---|
| 39 | id_=identifier; |
---|
| 40 | if (numberMembers_) { |
---|
| 41 | weights_ = new double[numberMembers_]; |
---|
[640] | 42 | which_ = new int[numberMembers_*numberLinks_]; |
---|
[121] | 43 | if (weights) { |
---|
| 44 | memcpy(weights_,weights,numberMembers_*sizeof(double)); |
---|
| 45 | } else { |
---|
| 46 | for (int i=0;i<numberMembers_;i++) |
---|
| 47 | weights_[i]=i; |
---|
| 48 | } |
---|
| 49 | // weights must be increasing |
---|
| 50 | int i; |
---|
[1902] | 51 | for (i=0;i<numberMembers_;i++) |
---|
| 52 | assert (i == 0 || weights_[i]>weights_[i-1]+1.0e-12); |
---|
[640] | 53 | for (i=0;i<numberMembers_*numberLinks_;i++) { |
---|
| 54 | which_[i]=first+i; |
---|
| 55 | } |
---|
[121] | 56 | } else { |
---|
| 57 | weights_ = NULL; |
---|
| 58 | } |
---|
| 59 | } |
---|
| 60 | |
---|
[640] | 61 | // Useful constructor (which are indices) |
---|
| 62 | CbcLink::CbcLink (CbcModel * model, int numberMembers, |
---|
| 63 | int numberLinks, int sosType, const int * which , const double * weights, int identifier) |
---|
| 64 | : CbcObject(model), |
---|
| 65 | numberMembers_(numberMembers), |
---|
| 66 | numberLinks_(numberLinks), |
---|
| 67 | which_(NULL), |
---|
| 68 | sosType_(sosType) |
---|
| 69 | { |
---|
| 70 | id_=identifier; |
---|
| 71 | if (numberMembers_) { |
---|
| 72 | weights_ = new double[numberMembers_]; |
---|
| 73 | which_ = new int[numberMembers_*numberLinks_]; |
---|
| 74 | if (weights) { |
---|
| 75 | memcpy(weights_,weights,numberMembers_*sizeof(double)); |
---|
| 76 | } else { |
---|
| 77 | for (int i=0;i<numberMembers_;i++) |
---|
| 78 | weights_[i]=i; |
---|
| 79 | } |
---|
| 80 | // weights must be increasing |
---|
| 81 | int i; |
---|
[1902] | 82 | for (i=0;i<numberMembers_;i++) |
---|
| 83 | assert (i == 0 || weights_[i]>weights_[i-1]+1.0e-12); |
---|
[640] | 84 | for (i=0;i<numberMembers_*numberLinks_;i++) { |
---|
| 85 | which_[i]= which[i]; |
---|
| 86 | } |
---|
| 87 | } else { |
---|
| 88 | weights_ = NULL; |
---|
| 89 | } |
---|
| 90 | } |
---|
| 91 | |
---|
[121] | 92 | // Copy constructor |
---|
| 93 | CbcLink::CbcLink ( const CbcLink & rhs) |
---|
| 94 | :CbcObject(rhs) |
---|
| 95 | { |
---|
| 96 | numberMembers_ = rhs.numberMembers_; |
---|
| 97 | numberLinks_ = rhs.numberLinks_; |
---|
[640] | 98 | sosType_ = rhs.sosType_; |
---|
[121] | 99 | if (numberMembers_) { |
---|
[640] | 100 | weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_); |
---|
| 101 | which_ = CoinCopyOfArray(rhs.which_,numberMembers_*numberLinks_); |
---|
[121] | 102 | } else { |
---|
| 103 | weights_ = NULL; |
---|
[640] | 104 | which_ = NULL; |
---|
[121] | 105 | } |
---|
| 106 | } |
---|
| 107 | |
---|
| 108 | // Clone |
---|
| 109 | CbcObject * |
---|
| 110 | CbcLink::clone() const |
---|
| 111 | { |
---|
| 112 | return new CbcLink(*this); |
---|
| 113 | } |
---|
| 114 | |
---|
| 115 | // Assignment operator |
---|
| 116 | CbcLink & |
---|
| 117 | CbcLink::operator=( const CbcLink& rhs) |
---|
| 118 | { |
---|
| 119 | if (this!=&rhs) { |
---|
| 120 | CbcObject::operator=(rhs); |
---|
| 121 | delete [] weights_; |
---|
[640] | 122 | delete [] which_; |
---|
[121] | 123 | numberMembers_ = rhs.numberMembers_; |
---|
| 124 | numberLinks_ = rhs.numberLinks_; |
---|
[640] | 125 | sosType_ = rhs.sosType_; |
---|
[121] | 126 | if (numberMembers_) { |
---|
[640] | 127 | weights_ = CoinCopyOfArray(rhs.weights_,numberMembers_); |
---|
| 128 | which_ = CoinCopyOfArray(rhs.which_,numberMembers_*numberLinks_); |
---|
[121] | 129 | } else { |
---|
| 130 | weights_ = NULL; |
---|
[640] | 131 | which_ = NULL; |
---|
[121] | 132 | } |
---|
| 133 | } |
---|
| 134 | return *this; |
---|
| 135 | } |
---|
| 136 | |
---|
| 137 | // Destructor |
---|
| 138 | CbcLink::~CbcLink () |
---|
| 139 | { |
---|
| 140 | delete [] weights_; |
---|
[640] | 141 | delete [] which_; |
---|
[121] | 142 | } |
---|
| 143 | |
---|
| 144 | // Infeasibility - large is 0.5 |
---|
| 145 | double |
---|
| 146 | CbcLink::infeasibility(int & preferredWay) const |
---|
| 147 | { |
---|
| 148 | int j; |
---|
| 149 | int firstNonZero=-1; |
---|
| 150 | int lastNonZero = -1; |
---|
| 151 | OsiSolverInterface * solver = model_->solver(); |
---|
[134] | 152 | const double * solution = model_->testSolution(); |
---|
[640] | 153 | //const double * lower = solver->getColLower(); |
---|
[121] | 154 | const double * upper = solver->getColUpper(); |
---|
| 155 | double integerTolerance = |
---|
| 156 | model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
| 157 | double weight = 0.0; |
---|
| 158 | double sum =0.0; |
---|
| 159 | |
---|
| 160 | // check bounds etc |
---|
| 161 | double lastWeight=-1.0e100; |
---|
[640] | 162 | int base=0; |
---|
[121] | 163 | for (j=0;j<numberMembers_;j++) { |
---|
| 164 | for (int k=0;k<numberLinks_;k++) { |
---|
[640] | 165 | int iColumn = which_[base+k]; |
---|
| 166 | //if (lower[iColumn]) |
---|
| 167 | //throw CoinError("Non zero lower bound in CBCLink","infeasibility","CbcLink"); |
---|
[121] | 168 | if (lastWeight>=weights_[j]-1.0e-7) |
---|
| 169 | throw CoinError("Weights too close together in CBCLink","infeasibility","CbcLink"); |
---|
| 170 | double value = CoinMax(0.0,solution[iColumn]); |
---|
| 171 | sum += value; |
---|
| 172 | if (value>integerTolerance&&upper[iColumn]) { |
---|
| 173 | // Possibly due to scaling a fixed variable might slip through |
---|
| 174 | if (value>upper[iColumn]+1.0e-8) { |
---|
| 175 | // Could change to #ifdef CBC_DEBUG |
---|
| 176 | #ifndef NDEBUG |
---|
| 177 | if (model_->messageHandler()->logLevel()>1) |
---|
| 178 | printf("** Variable %d (%d) has value %g and upper bound of %g\n", |
---|
| 179 | iColumn,j,value,upper[iColumn]); |
---|
| 180 | #endif |
---|
| 181 | } |
---|
[640] | 182 | value = CoinMin(value,upper[iColumn]); |
---|
| 183 | weight += weights_[j]*value; |
---|
[121] | 184 | if (firstNonZero<0) |
---|
| 185 | firstNonZero=j; |
---|
| 186 | lastNonZero=j; |
---|
| 187 | } |
---|
| 188 | } |
---|
| 189 | base += numberLinks_; |
---|
| 190 | } |
---|
[640] | 191 | double valueInfeasibility; |
---|
[121] | 192 | preferredWay=1; |
---|
[640] | 193 | if (lastNonZero-firstNonZero>=sosType_) { |
---|
[121] | 194 | // find where to branch |
---|
| 195 | assert (sum>0.0); |
---|
| 196 | weight /= sum; |
---|
[640] | 197 | valueInfeasibility = lastNonZero-firstNonZero+1; |
---|
| 198 | valueInfeasibility *= 0.5/((double) numberMembers_); |
---|
| 199 | //#define DISTANCE |
---|
| 200 | #ifdef DISTANCE |
---|
| 201 | assert (sosType_==1); // code up |
---|
| 202 | /* may still be satisfied. |
---|
| 203 | For LOS type 2 we might wish to move coding around |
---|
| 204 | and keep initial info in model_ for speed |
---|
| 205 | */ |
---|
| 206 | int iWhere; |
---|
| 207 | bool possible=false; |
---|
| 208 | for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) { |
---|
| 209 | if (fabs(weight-weights_[iWhere])<1.0e-8) { |
---|
| 210 | possible=true; |
---|
| 211 | break; |
---|
| 212 | } |
---|
| 213 | } |
---|
| 214 | if (possible) { |
---|
| 215 | // One could move some of this (+ arrays) into model_ |
---|
| 216 | const CoinPackedMatrix * matrix = solver->getMatrixByCol(); |
---|
| 217 | const double * element = matrix->getMutableElements(); |
---|
| 218 | const int * row = matrix->getIndices(); |
---|
| 219 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
---|
| 220 | const int * columnLength = matrix->getVectorLengths(); |
---|
| 221 | const double * rowSolution = solver->getRowActivity(); |
---|
| 222 | const double * rowLower = solver->getRowLower(); |
---|
| 223 | const double * rowUpper = solver->getRowUpper(); |
---|
| 224 | int numberRows = matrix->getNumRows(); |
---|
| 225 | double * array = new double [numberRows]; |
---|
| 226 | CoinZeroN(array,numberRows); |
---|
| 227 | int * which = new int [numberRows]; |
---|
| 228 | int n=0; |
---|
| 229 | int base=numberLinks_*firstNonZero; |
---|
| 230 | for (j=firstNonZero;j<=lastNonZero;j++) { |
---|
| 231 | for (int k=0;k<numberLinks_;k++) { |
---|
| 232 | int iColumn = which_[base+k]; |
---|
| 233 | double value = CoinMax(0.0,solution[iColumn]); |
---|
| 234 | if (value>integerTolerance&&upper[iColumn]) { |
---|
| 235 | value = CoinMin(value,upper[iColumn]); |
---|
| 236 | for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { |
---|
| 237 | int iRow = row[j]; |
---|
| 238 | double a = array[iRow]; |
---|
| 239 | if (a) { |
---|
| 240 | a += value*element[j]; |
---|
| 241 | if (!a) |
---|
| 242 | a = 1.0e-100; |
---|
| 243 | } else { |
---|
| 244 | which[n++]=iRow; |
---|
| 245 | a=value*element[j]; |
---|
| 246 | assert (a); |
---|
| 247 | } |
---|
| 248 | array[iRow]=a; |
---|
| 249 | } |
---|
| 250 | } |
---|
| 251 | } |
---|
| 252 | base += numberLinks_; |
---|
| 253 | } |
---|
| 254 | base=numberLinks_*iWhere; |
---|
| 255 | for (int k=0;k<numberLinks_;k++) { |
---|
| 256 | int iColumn = which_[base+k]; |
---|
| 257 | const double value = 1.0; |
---|
| 258 | for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { |
---|
| 259 | int iRow = row[j]; |
---|
| 260 | double a = array[iRow]; |
---|
| 261 | if (a) { |
---|
| 262 | a -= value*element[j]; |
---|
| 263 | if (!a) |
---|
| 264 | a = 1.0e-100; |
---|
| 265 | } else { |
---|
| 266 | which[n++]=iRow; |
---|
| 267 | a=-value*element[j]; |
---|
| 268 | assert (a); |
---|
| 269 | } |
---|
| 270 | array[iRow]=a; |
---|
| 271 | } |
---|
| 272 | } |
---|
| 273 | for (j=0;j<n;j++) { |
---|
| 274 | int iRow = which[j]; |
---|
| 275 | // moving to point will increase row solution by this |
---|
| 276 | double distance = array[iRow]; |
---|
| 277 | if (distance>1.0e-8) { |
---|
| 278 | if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) { |
---|
| 279 | possible=false; |
---|
| 280 | break; |
---|
| 281 | } |
---|
| 282 | } else if (distance<-1.0e-8) { |
---|
| 283 | if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) { |
---|
| 284 | possible=false; |
---|
| 285 | break; |
---|
| 286 | } |
---|
| 287 | } |
---|
| 288 | } |
---|
| 289 | for (j=0;j<n;j++) |
---|
| 290 | array[which[j]]=0.0; |
---|
| 291 | delete [] array; |
---|
| 292 | delete [] which; |
---|
| 293 | if (possible) { |
---|
| 294 | valueInfeasibility=0.0; |
---|
| 295 | printf("possible %d %d %d\n",firstNonZero,lastNonZero,iWhere); |
---|
| 296 | } |
---|
| 297 | } |
---|
| 298 | #endif |
---|
[121] | 299 | } else { |
---|
[640] | 300 | valueInfeasibility = 0.0; // satisfied |
---|
[121] | 301 | } |
---|
[640] | 302 | return valueInfeasibility; |
---|
[121] | 303 | } |
---|
| 304 | |
---|
| 305 | // This looks at solution and sets bounds to contain solution |
---|
| 306 | void |
---|
| 307 | CbcLink::feasibleRegion() |
---|
| 308 | { |
---|
| 309 | int j; |
---|
| 310 | int firstNonZero=-1; |
---|
| 311 | int lastNonZero = -1; |
---|
| 312 | OsiSolverInterface * solver = model_->solver(); |
---|
[134] | 313 | const double * solution = model_->testSolution(); |
---|
[121] | 314 | const double * upper = solver->getColUpper(); |
---|
| 315 | double integerTolerance = |
---|
| 316 | model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
| 317 | double weight = 0.0; |
---|
| 318 | double sum =0.0; |
---|
| 319 | |
---|
[640] | 320 | int base=0; |
---|
[121] | 321 | for (j=0;j<numberMembers_;j++) { |
---|
| 322 | for (int k=0;k<numberLinks_;k++) { |
---|
[640] | 323 | int iColumn = which_[base+k]; |
---|
[121] | 324 | double value = CoinMax(0.0,solution[iColumn]); |
---|
| 325 | sum += value; |
---|
| 326 | if (value>integerTolerance&&upper[iColumn]) { |
---|
| 327 | weight += weights_[j]*value; |
---|
| 328 | if (firstNonZero<0) |
---|
| 329 | firstNonZero=j; |
---|
| 330 | lastNonZero=j; |
---|
| 331 | } |
---|
| 332 | } |
---|
| 333 | base += numberLinks_; |
---|
| 334 | } |
---|
[640] | 335 | #ifdef DISTANCE |
---|
| 336 | if (lastNonZero-firstNonZero>sosType_-1) { |
---|
| 337 | /* may still be satisfied. |
---|
| 338 | For LOS type 2 we might wish to move coding around |
---|
| 339 | and keep initial info in model_ for speed |
---|
| 340 | */ |
---|
| 341 | int iWhere; |
---|
| 342 | bool possible=false; |
---|
| 343 | for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) { |
---|
| 344 | if (fabs(weight-weights_[iWhere])<1.0e-8) { |
---|
| 345 | possible=true; |
---|
| 346 | break; |
---|
| 347 | } |
---|
| 348 | } |
---|
| 349 | if (possible) { |
---|
| 350 | // One could move some of this (+ arrays) into model_ |
---|
| 351 | const CoinPackedMatrix * matrix = solver->getMatrixByCol(); |
---|
| 352 | const double * element = matrix->getMutableElements(); |
---|
| 353 | const int * row = matrix->getIndices(); |
---|
| 354 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
---|
| 355 | const int * columnLength = matrix->getVectorLengths(); |
---|
| 356 | const double * rowSolution = solver->getRowActivity(); |
---|
| 357 | const double * rowLower = solver->getRowLower(); |
---|
| 358 | const double * rowUpper = solver->getRowUpper(); |
---|
| 359 | int numberRows = matrix->getNumRows(); |
---|
| 360 | double * array = new double [numberRows]; |
---|
| 361 | CoinZeroN(array,numberRows); |
---|
| 362 | int * which = new int [numberRows]; |
---|
| 363 | int n=0; |
---|
| 364 | int base=numberLinks_*firstNonZero; |
---|
| 365 | for (j=firstNonZero;j<=lastNonZero;j++) { |
---|
| 366 | for (int k=0;k<numberLinks_;k++) { |
---|
| 367 | int iColumn = which_[base+k]; |
---|
| 368 | double value = CoinMax(0.0,solution[iColumn]); |
---|
| 369 | if (value>integerTolerance&&upper[iColumn]) { |
---|
| 370 | value = CoinMin(value,upper[iColumn]); |
---|
| 371 | for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { |
---|
| 372 | int iRow = row[j]; |
---|
| 373 | double a = array[iRow]; |
---|
| 374 | if (a) { |
---|
| 375 | a += value*element[j]; |
---|
| 376 | if (!a) |
---|
| 377 | a = 1.0e-100; |
---|
| 378 | } else { |
---|
| 379 | which[n++]=iRow; |
---|
| 380 | a=value*element[j]; |
---|
| 381 | assert (a); |
---|
| 382 | } |
---|
| 383 | array[iRow]=a; |
---|
| 384 | } |
---|
| 385 | } |
---|
| 386 | } |
---|
| 387 | base += numberLinks_; |
---|
| 388 | } |
---|
| 389 | base=numberLinks_*iWhere; |
---|
| 390 | for (int k=0;k<numberLinks_;k++) { |
---|
| 391 | int iColumn = which_[base+k]; |
---|
| 392 | const double value = 1.0; |
---|
| 393 | for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) { |
---|
| 394 | int iRow = row[j]; |
---|
| 395 | double a = array[iRow]; |
---|
| 396 | if (a) { |
---|
| 397 | a -= value*element[j]; |
---|
| 398 | if (!a) |
---|
| 399 | a = 1.0e-100; |
---|
| 400 | } else { |
---|
| 401 | which[n++]=iRow; |
---|
| 402 | a=-value*element[j]; |
---|
| 403 | assert (a); |
---|
| 404 | } |
---|
| 405 | array[iRow]=a; |
---|
| 406 | } |
---|
| 407 | } |
---|
| 408 | for (j=0;j<n;j++) { |
---|
| 409 | int iRow = which[j]; |
---|
| 410 | // moving to point will increase row solution by this |
---|
| 411 | double distance = array[iRow]; |
---|
| 412 | if (distance>1.0e-8) { |
---|
| 413 | if (distance+rowSolution[iRow]>rowUpper[iRow]+1.0e-8) { |
---|
| 414 | possible=false; |
---|
| 415 | break; |
---|
| 416 | } |
---|
| 417 | } else if (distance<-1.0e-8) { |
---|
| 418 | if (distance+rowSolution[iRow]<rowLower[iRow]-1.0e-8) { |
---|
| 419 | possible=false; |
---|
| 420 | break; |
---|
| 421 | } |
---|
| 422 | } |
---|
| 423 | } |
---|
| 424 | for (j=0;j<n;j++) |
---|
| 425 | array[which[j]]=0.0; |
---|
| 426 | delete [] array; |
---|
| 427 | delete [] which; |
---|
| 428 | if (possible) { |
---|
| 429 | printf("possible feas region %d %d %d\n",firstNonZero,lastNonZero,iWhere); |
---|
| 430 | firstNonZero=iWhere; |
---|
| 431 | lastNonZero=iWhere; |
---|
| 432 | } |
---|
| 433 | } |
---|
| 434 | } |
---|
| 435 | #else |
---|
| 436 | assert (lastNonZero-firstNonZero<sosType_) ; |
---|
| 437 | #endif |
---|
| 438 | base=0; |
---|
[121] | 439 | for (j=0;j<firstNonZero;j++) { |
---|
| 440 | for (int k=0;k<numberLinks_;k++) { |
---|
[640] | 441 | int iColumn = which_[base+k]; |
---|
[121] | 442 | solver->setColUpper(iColumn,0.0); |
---|
| 443 | } |
---|
| 444 | base += numberLinks_; |
---|
| 445 | } |
---|
| 446 | // skip |
---|
| 447 | base += numberLinks_; |
---|
| 448 | for (j=lastNonZero+1;j<numberMembers_;j++) { |
---|
| 449 | for (int k=0;k<numberLinks_;k++) { |
---|
[640] | 450 | int iColumn = which_[base+k]; |
---|
[121] | 451 | solver->setColUpper(iColumn,0.0); |
---|
| 452 | } |
---|
| 453 | base += numberLinks_; |
---|
| 454 | } |
---|
| 455 | } |
---|
| 456 | |
---|
| 457 | |
---|
| 458 | // Creates a branching object |
---|
| 459 | CbcBranchingObject * |
---|
[1902] | 460 | CbcLink::createCbcBranch(OsiSolverInterface * /*solver*/, const OsiBranchingInformation * /*info*/, int way) |
---|
[121] | 461 | { |
---|
| 462 | int j; |
---|
[134] | 463 | const double * solution = model_->testSolution(); |
---|
[121] | 464 | double integerTolerance = |
---|
| 465 | model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
| 466 | OsiSolverInterface * solver = model_->solver(); |
---|
| 467 | const double * upper = solver->getColUpper(); |
---|
| 468 | int firstNonFixed=-1; |
---|
| 469 | int lastNonFixed=-1; |
---|
| 470 | int firstNonZero=-1; |
---|
| 471 | int lastNonZero = -1; |
---|
| 472 | double weight = 0.0; |
---|
| 473 | double sum =0.0; |
---|
[640] | 474 | int base=0; |
---|
[121] | 475 | for (j=0;j<numberMembers_;j++) { |
---|
| 476 | for (int k=0;k<numberLinks_;k++) { |
---|
[640] | 477 | int iColumn = which_[base+k]; |
---|
[121] | 478 | if (upper[iColumn]) { |
---|
| 479 | double value = CoinMax(0.0,solution[iColumn]); |
---|
| 480 | sum += value; |
---|
| 481 | if (firstNonFixed<0) |
---|
| 482 | firstNonFixed=j; |
---|
| 483 | lastNonFixed=j; |
---|
| 484 | if (value>integerTolerance) { |
---|
| 485 | weight += weights_[j]*value; |
---|
| 486 | if (firstNonZero<0) |
---|
| 487 | firstNonZero=j; |
---|
| 488 | lastNonZero=j; |
---|
| 489 | } |
---|
| 490 | } |
---|
| 491 | } |
---|
| 492 | base += numberLinks_; |
---|
| 493 | } |
---|
[640] | 494 | assert (lastNonZero-firstNonZero>=sosType_) ; |
---|
[121] | 495 | // find where to branch |
---|
| 496 | assert (sum>0.0); |
---|
| 497 | weight /= sum; |
---|
| 498 | int iWhere; |
---|
| 499 | double separator=0.0; |
---|
| 500 | for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) |
---|
| 501 | if (weight<weights_[iWhere+1]) |
---|
| 502 | break; |
---|
[640] | 503 | if (sosType_==1) { |
---|
| 504 | // SOS 1 |
---|
| 505 | separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]); |
---|
| 506 | } else { |
---|
| 507 | // SOS 2 |
---|
| 508 | if (iWhere==firstNonFixed) |
---|
| 509 | iWhere++;; |
---|
| 510 | if (iWhere==lastNonFixed-1) |
---|
| 511 | iWhere = lastNonFixed-2; |
---|
| 512 | separator = weights_[iWhere+1]; |
---|
| 513 | } |
---|
[121] | 514 | // create object |
---|
| 515 | CbcBranchingObject * branch; |
---|
| 516 | branch = new CbcLinkBranchingObject(model_,this,way,separator); |
---|
| 517 | return branch; |
---|
| 518 | } |
---|
| 519 | // Useful constructor |
---|
| 520 | CbcLinkBranchingObject::CbcLinkBranchingObject (CbcModel * model, |
---|
| 521 | const CbcLink * set, |
---|
| 522 | int way , |
---|
| 523 | double separator) |
---|
| 524 | :CbcBranchingObject(model,set->id(),way,0.5) |
---|
| 525 | { |
---|
| 526 | set_ = set; |
---|
| 527 | separator_ = separator; |
---|
| 528 | } |
---|
| 529 | |
---|
| 530 | // Copy constructor |
---|
| 531 | CbcLinkBranchingObject::CbcLinkBranchingObject ( const CbcLinkBranchingObject & rhs) :CbcBranchingObject(rhs) |
---|
| 532 | { |
---|
| 533 | set_=rhs.set_; |
---|
| 534 | separator_ = rhs.separator_; |
---|
| 535 | } |
---|
| 536 | |
---|
| 537 | // Assignment operator |
---|
| 538 | CbcLinkBranchingObject & |
---|
| 539 | CbcLinkBranchingObject::operator=( const CbcLinkBranchingObject& rhs) |
---|
| 540 | { |
---|
| 541 | if (this != &rhs) { |
---|
| 542 | CbcBranchingObject::operator=(rhs); |
---|
| 543 | set_=rhs.set_; |
---|
| 544 | separator_ = rhs.separator_; |
---|
| 545 | } |
---|
| 546 | return *this; |
---|
| 547 | } |
---|
| 548 | CbcBranchingObject * |
---|
| 549 | CbcLinkBranchingObject::clone() const |
---|
| 550 | { |
---|
| 551 | return (new CbcLinkBranchingObject(*this)); |
---|
| 552 | } |
---|
| 553 | |
---|
| 554 | |
---|
| 555 | // Destructor |
---|
| 556 | CbcLinkBranchingObject::~CbcLinkBranchingObject () |
---|
| 557 | { |
---|
| 558 | } |
---|
| 559 | double |
---|
[640] | 560 | CbcLinkBranchingObject::branch() |
---|
[121] | 561 | { |
---|
[640] | 562 | decrementNumberBranchesLeft(); |
---|
[121] | 563 | int numberMembers = set_->numberMembers(); |
---|
| 564 | int numberLinks = set_->numberLinks(); |
---|
| 565 | const double * weights = set_->weights(); |
---|
[640] | 566 | const int * which = set_->which(); |
---|
[121] | 567 | OsiSolverInterface * solver = model_->solver(); |
---|
| 568 | //const double * lower = solver->getColLower(); |
---|
| 569 | //const double * upper = solver->getColUpper(); |
---|
| 570 | // *** for way - up means fix all those in down section |
---|
| 571 | if (way_<0) { |
---|
| 572 | int i; |
---|
| 573 | for ( i=0;i<numberMembers;i++) { |
---|
| 574 | if (weights[i] > separator_) |
---|
| 575 | break; |
---|
| 576 | } |
---|
| 577 | assert (i<numberMembers); |
---|
[640] | 578 | int base=i*numberLinks;; |
---|
[121] | 579 | for (;i<numberMembers;i++) { |
---|
| 580 | for (int k=0;k<numberLinks;k++) { |
---|
[640] | 581 | int iColumn = which[base+k]; |
---|
[121] | 582 | solver->setColUpper(iColumn,0.0); |
---|
| 583 | } |
---|
| 584 | base += numberLinks; |
---|
| 585 | } |
---|
| 586 | way_=1; // Swap direction |
---|
| 587 | } else { |
---|
| 588 | int i; |
---|
[640] | 589 | int base=0; |
---|
[121] | 590 | for ( i=0;i<numberMembers;i++) { |
---|
| 591 | if (weights[i] >= separator_) { |
---|
| 592 | break; |
---|
| 593 | } else { |
---|
| 594 | for (int k=0;k<numberLinks;k++) { |
---|
[640] | 595 | int iColumn = which[base+k]; |
---|
[121] | 596 | solver->setColUpper(iColumn,0.0); |
---|
| 597 | } |
---|
| 598 | base += numberLinks; |
---|
| 599 | } |
---|
| 600 | } |
---|
| 601 | assert (i<numberMembers); |
---|
| 602 | way_=-1; // Swap direction |
---|
| 603 | } |
---|
| 604 | return 0.0; |
---|
| 605 | } |
---|
| 606 | // Print what would happen |
---|
| 607 | void |
---|
[640] | 608 | CbcLinkBranchingObject::print() |
---|
[121] | 609 | { |
---|
| 610 | int numberMembers = set_->numberMembers(); |
---|
| 611 | int numberLinks = set_->numberLinks(); |
---|
| 612 | const double * weights = set_->weights(); |
---|
[640] | 613 | const int * which = set_->which(); |
---|
[121] | 614 | OsiSolverInterface * solver = model_->solver(); |
---|
| 615 | const double * upper = solver->getColUpper(); |
---|
| 616 | int first=numberMembers; |
---|
| 617 | int last=-1; |
---|
| 618 | int numberFixed=0; |
---|
| 619 | int numberOther=0; |
---|
| 620 | int i; |
---|
[640] | 621 | int base=0; |
---|
[121] | 622 | for ( i=0;i<numberMembers;i++) { |
---|
| 623 | for (int k=0;k<numberLinks;k++) { |
---|
[640] | 624 | int iColumn = which[base+k]; |
---|
[121] | 625 | double bound = upper[iColumn]; |
---|
| 626 | if (bound) { |
---|
| 627 | first = CoinMin(first,i); |
---|
| 628 | last = CoinMax(last,i); |
---|
| 629 | } |
---|
| 630 | } |
---|
| 631 | base += numberLinks; |
---|
| 632 | } |
---|
| 633 | // *** for way - up means fix all those in down section |
---|
[640] | 634 | base=0; |
---|
[121] | 635 | if (way_<0) { |
---|
| 636 | printf("SOS Down"); |
---|
| 637 | for ( i=0;i<numberMembers;i++) { |
---|
[640] | 638 | if (weights[i] > separator_) |
---|
| 639 | break; |
---|
[121] | 640 | for (int k=0;k<numberLinks;k++) { |
---|
[640] | 641 | int iColumn = which[base+k]; |
---|
[121] | 642 | double bound = upper[iColumn]; |
---|
[640] | 643 | if (bound) |
---|
[121] | 644 | numberOther++; |
---|
| 645 | } |
---|
| 646 | base += numberLinks; |
---|
| 647 | } |
---|
| 648 | assert (i<numberMembers); |
---|
| 649 | for (;i<numberMembers;i++) { |
---|
| 650 | for (int k=0;k<numberLinks;k++) { |
---|
[640] | 651 | int iColumn = which[base+k]; |
---|
[121] | 652 | double bound = upper[iColumn]; |
---|
| 653 | if (bound) |
---|
| 654 | numberFixed++; |
---|
| 655 | } |
---|
| 656 | base += numberLinks; |
---|
| 657 | } |
---|
| 658 | } else { |
---|
| 659 | printf("SOS Up"); |
---|
| 660 | for ( i=0;i<numberMembers;i++) { |
---|
[640] | 661 | if (weights[i] >= separator_) |
---|
| 662 | break; |
---|
[121] | 663 | for (int k=0;k<numberLinks;k++) { |
---|
[640] | 664 | int iColumn = which[base+k]; |
---|
[121] | 665 | double bound = upper[iColumn]; |
---|
[640] | 666 | if (bound) |
---|
[121] | 667 | numberFixed++; |
---|
| 668 | } |
---|
| 669 | base += numberLinks; |
---|
| 670 | } |
---|
| 671 | assert (i<numberMembers); |
---|
| 672 | for (;i<numberMembers;i++) { |
---|
| 673 | for (int k=0;k<numberLinks;k++) { |
---|
[640] | 674 | int iColumn = which[base+k]; |
---|
[121] | 675 | double bound = upper[iColumn]; |
---|
| 676 | if (bound) |
---|
| 677 | numberOther++; |
---|
| 678 | } |
---|
| 679 | base += numberLinks; |
---|
| 680 | } |
---|
| 681 | } |
---|
| 682 | assert ((numberFixed%numberLinks)==0); |
---|
[640] | 683 | assert ((numberOther%numberLinks)==0); |
---|
[121] | 684 | printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n", |
---|
| 685 | separator_,first,weights[first],last,weights[last],numberFixed/numberLinks, |
---|
| 686 | numberOther/numberLinks); |
---|
| 687 | } |
---|
[932] | 688 | /** Compare the \c this with \c brObj. \c this and \c brObj must be os the |
---|
| 689 | same type and must have the same original object, but they may have |
---|
| 690 | different feasible regions. |
---|
| 691 | Return the appropriate CbcRangeCompare value (first argument being the |
---|
| 692 | sub/superset if that's the case). In case of overlap (and if \c |
---|
| 693 | replaceIfOverlap is true) replace the current branching object with one |
---|
| 694 | whose feasible region is the overlap. |
---|
| 695 | */ |
---|
| 696 | CbcRangeCompare |
---|
| 697 | CbcLinkBranchingObject::compareBranchingObject |
---|
| 698 | (const CbcBranchingObject* brObj, const bool replaceIfOverlap) |
---|
| 699 | { |
---|
| 700 | throw("must implement"); |
---|
| 701 | } |
---|