1  // Copyright (C) 2002, International Business Machines 

2  // Corporation and others. All Rights Reserved. 

3  #if defined(_MSC_VER) 

4  // Turn off compiler warning about long names 

5  # pragma warning(disable:4786) 

6  #endif 

7  #include <string> 

8  //#define CBC_DEBUG 1 

9  //#define CHECK_CUT_COUNTS 

10  //#define CHECK_NODE_FULL 

11  //#define NODE_LOG 

12  #include <cassert> 

13  #include <cmath> 

14  #include <cfloat> 

15  

16  #ifdef CBC_USE_CLP 

17  // include Presolve from Clp 

18  #include "ClpPresolve.hpp" 

19  #include "OsiClpSolverInterface.hpp" 

20  #endif 

21  

22  #include "CbcEventHandler.hpp" 

23  

24  #include "OsiSolverInterface.hpp" 

25  #include "OsiAuxInfo.hpp" 

26  #include "OsiSolverBranch.hpp" 

27  #include "CoinWarmStartBasis.hpp" 

28  #include "CoinPackedMatrix.hpp" 

29  #include "CoinHelperFunctions.hpp" 

30  #include "CbcBranchActual.hpp" 

31  #include "CbcBranchDynamic.hpp" 

32  #include "CbcHeuristic.hpp" 

33  #include "CbcModel.hpp" 

34  #include "CbcTreeLocal.hpp" 

35  #include "CbcStatistics.hpp" 

36  #include "CbcStrategy.hpp" 

37  #include "CbcMessage.hpp" 

38  #include "OsiRowCut.hpp" 

39  #include "OsiColCut.hpp" 

40  #include "OsiRowCutDebugger.hpp" 

41  #include "OsiCuts.hpp" 

42  #include "CbcCountRowCut.hpp" 

43  #include "CbcCutGenerator.hpp" 

44  #include "CbcFeasibilityBase.hpp" 

45  // include Probing 

46  #include "CglProbing.hpp" 

47  // include preprocessing 

48  #include "CglPreProcess.hpp" 

49  

50  #include "CoinTime.hpp" 

51  

52  #include "CbcCompareActual.hpp" 

53  #include "CbcTree.hpp" 

54  /* Various functions local to CbcModel.cpp */ 

55  

56  namespace { 

57  

58  // 

59  // Returns the greatest common denominator of two 

60  // positive integers, a and b, found using Euclid's algorithm 

61  // 

62  static int gcd(int a, int b) 

63  { 

64  int remainder = 1; 

65  // make sure a<=b (will always remain so) 

66  if(a > b) { 

67  // Swap a and b 

68  int temp = a; 

69  a = b; 

70  b = temp; 

71  } 

72  // if zero then gcd is nonzero (zero may occur in rhs of packed) 

73  if (!a) { 

74  if (b) { 

75  return b; 

76  } else { 

77  printf("**** gcd given two zeros!!\n"); 

78  abort(); 

79  } 

80  } 

81  while (remainder) { 

82  remainder = b % a; 

83  b = a; 

84  a = remainder; 

85  } 

86  return b; 

87  } 

88  

89  

90  

91  #ifdef CHECK_NODE_FULL 

92  

93  /* 

94  Routine to verify that tree linkage is correct. The invariant that is tested 

95  is 

96  

97  reference count = (number of actual references) + (number of branches left) 

98  

99  The routine builds a set of paired arrays, info and count, by traversing the 

100  tree. Each CbcNodeInfo is recorded in info, and the number of times it is 

101  referenced (via the parent field) is recorded in count. Then a final check is 

102  made to see if the numberPointingToThis_ field agrees. 

103  */ 

104  

105  void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model) 

106  

107  { printf("*** CHECKING tree after %d nodes\n",model.getNodeCount()) ; 

108  

109  int j ; 

110  int nNodes = branchingTree>size() ; 

111  # define MAXINFO 1000 

112  int *count = new int [MAXINFO] ; 

113  CbcNodeInfo **info = new CbcNodeInfo*[MAXINFO] ; 

114  int nInfo = 0 ; 

115  /* 

116  Collect all CbcNodeInfo objects in info, by starting from each live node and 

117  traversing back to the root. Nodes in the live set should have unexplored 

118  branches remaining. 

119  

120  TODO: The `while (nodeInfo)' loop could be made to break on reaching a 

121  common ancester (nodeInfo is found in info[k]). Alternatively, the 

122  check could change to signal an error if nodeInfo is not found above a 

123  common ancestor. 

124  */ 

125  for (j = 0 ; j < nNodes ; j++) 

126  { CbcNode *node = branchingTree>nodePointer(j) ; 

127  if (!node) 

128  continue; 

129  CbcNodeInfo *nodeInfo = node>nodeInfo() ; 

130  int change = node>nodeInfo()>numberBranchesLeft() ; 

131  assert(change) ; 

132  while (nodeInfo) 

133  { int k ; 

134  for (k = 0 ; k < nInfo ; k++) 

135  { if (nodeInfo == info[k]) break ; } 

136  if (k == nInfo) 

137  { assert(nInfo < MAXINFO) ; 

138  nInfo++ ; 

139  info[k] = nodeInfo ; 

140  count[k] = 0 ; } 

141  nodeInfo = nodeInfo>parent() ; } } 

142  /* 

143  Walk the info array. For each nodeInfo, look up its parent in info and 

144  increment the corresponding count. 

145  */ 

146  for (j = 0 ; j < nInfo ; j++) 

147  { CbcNodeInfo *nodeInfo = info[j] ; 

148  nodeInfo = nodeInfo>parent() ; 

149  if (nodeInfo) 

150  { int k ; 

151  for (k = 0 ; k < nInfo ; k++) 

152  { if (nodeInfo == info[k]) break ; } 

153  assert (k < nInfo) ; 

154  count[k]++ ; } } 

155  /* 

156  Walk the info array one more time and check that the invariant holds. The 

157  number of references (numberPointingToThis()) should equal the sum of the 

158  number of actual references (held in count[]) plus the number of potential 

159  references (unexplored branches, numberBranchesLeft()). 

160  */ 

161  for (j = 0;j < nInfo;j++) { 

162  CbcNodeInfo * nodeInfo = info[j] ; 

163  if (nodeInfo) { 

164  int k ; 

165  for (k = 0;k < nInfo;k++) 

166  if (nodeInfo == info[k]) 

167  break ; 

168  printf("Nodeinfo %x  %d left, %d count\n", 

169  nodeInfo, 

170  nodeInfo>numberBranchesLeft(), 

171  nodeInfo>numberPointingToThis()) ; 

172  assert(nodeInfo>numberPointingToThis() == 

173  count[k]+nodeInfo>numberBranchesLeft()) ; } } 

174  

175  delete [] count ; 

176  delete [] info ; 

177  

178  return ; } 

179  

180  #endif /* CHECK_NODE_FULL */ 

181  

182  

183  

184  #ifdef CHECK_CUT_COUNTS 

185  

186  /* 

187  Routine to verify that cut reference counts are correct. 

188  */ 

189  void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model) 

190  

191  { printf("*** CHECKING cuts after %d nodes\n",model.getNodeCount()) ; 

192  

193  int j ; 

194  int nNodes = branchingTree>size() ; 

195  

196  /* 

197  cut.tempNumber_ exists for the purpose of doing this verification. Clear it 

198  in all cuts. We traverse the tree by starting from each live node and working 

199  back to the root. At each CbcNodeInfo, check for cuts. 

200  */ 

201  for (j = 0 ; j < nNodes ; j++) 

202  { CbcNode *node = branchingTree>nodePointer(j) ; 

203  CbcNodeInfo * nodeInfo = node>nodeInfo() ; 

204  assert (node>nodeInfo()>numberBranchesLeft()) ; 

205  while (nodeInfo) 

206  { int k ; 

207  for (k = 0 ; k < nodeInfo>numberCuts() ; k++) 

208  { CbcCountRowCut *cut = nodeInfo>cuts()[k] ; 

209  if (cut) cut>tempNumber_ = 0; } 

210  nodeInfo = nodeInfo>parent() ; } } 

211  /* 

212  Walk the live set again, this time collecting the list of cuts in use at each 

213  node. addCuts1 will collect the cuts in model.addedCuts_. Take into account 

214  that when we recreate the basis for a node, we compress out the slack cuts. 

215  */ 

216  for (j = 0 ; j < nNodes ; j++) 

217  { CoinWarmStartBasis *debugws = model.getEmptyBasis() ; 

218  CbcNode *node = branchingTree>nodePointer(j) ; 

219  CbcNodeInfo *nodeInfo = node>nodeInfo(); 

220  int change = node>nodeInfo()>numberBranchesLeft() ; 

221  printf("Node %d %x (info %x) var %d way %d obj %g",j,node, 

222  node>nodeInfo(),node>variable(),node>way(), 

223  node>objectiveValue()) ; 

224  

225  model.addCuts1(node,debugws) ; 

226  

227  int i ; 

228  int numberRowsAtContinuous = model.numberRowsAtContinuous() ; 

229  CbcCountRowCut **addedCuts = model.addedCuts() ; 

230  for (i = 0 ; i < model.currentNumberCuts() ; i++) 

231  { CoinWarmStartBasis::Status status = 

232  debugws>getArtifStatus(i+numberRowsAtContinuous) ; 

233  if (status != CoinWarmStartBasis::basic && addedCuts[i]) 

234  { addedCuts[i]>tempNumber_ += change ; } } 

235  

236  while (nodeInfo) 

237  { nodeInfo = nodeInfo>parent() ; 

238  if (nodeInfo) printf(" > %x",nodeInfo); } 

239  printf("\n") ; 

240  delete debugws ; } 

241  /* 

242  The moment of truth: We've tallied up the references by direct scan of the 

243  search tree. Check for agreement with the count in the cut. 

244  

245  TODO: Rewrite to check and print mismatch only when tempNumber_ == 0? 

246  */ 

247  for (j = 0 ; j < nNodes ; j++) 

248  { CbcNode *node = branchingTree>nodePointer(j) ; 

249  CbcNodeInfo *nodeInfo = node>nodeInfo(); 

250  while (nodeInfo) 

251  { int k ; 

252  for (k = 0 ; k < nodeInfo>numberCuts() ; k++) 

253  { CbcCountRowCut *cut = nodeInfo>cuts()[k] ; 

254  if (cut && cut>tempNumber_ >= 0) 

255  { if (cut>tempNumber_ != cut>numberPointingToThis()) 

256  printf("mismatch %x %d %x %d %d\n",nodeInfo,k, 

257  cut,cut>tempNumber_,cut>numberPointingToThis()) ; 

258  else 

259  printf(" match %x %d %x %d %d\n", nodeInfo,k, 

260  cut,cut>tempNumber_,cut>numberPointingToThis()) ; 

261  cut>tempNumber_ = 1 ; } } 

262  nodeInfo = nodeInfo>parent() ; } } 

263  

264  return ; } 

265  

266  #endif /* CHECK_CUT_COUNTS */ 

267  

268  } 

269  

270  /* End unnamed namespace for CbcModel.cpp */ 

271  

272  

273  

274  void 

275  CbcModel::analyzeObjective () 

276  /* 

277  Try to find a minimum change in the objective function. The first scan 

278  checks that there are no continuous variables with nonzero coefficients, 

279  and grabs the largest objective coefficient associated with an unfixed 

280  integer variable. The second scan attempts to scale up the objective 

281  coefficients to a point where they are sufficiently close to integer that 

282  we can pretend they are integer, and calculate a gcd over the coefficients 

283  of interest. This will be the minimum increment for the scaled coefficients. 

284  The final action is to scale the increment back for the original coefficients 

285  and install it, if it's better than the existing value. 

286  

287  John's note: We could do better than this. 

288  

289  John's second note  apologies for changing s to z 

290  */ 

291  { const double *objective = getObjCoefficients() ; 

292  const double *lower = getColLower() ; 

293  const double *upper = getColUpper() ; 

294  /* 

295  Take a first scan to see if there are unfixed continuous variables in the 

296  objective. If so, the minimum objective change could be arbitrarily small. 

297  Also pick off the maximum coefficient of an unfixed integer variable. 

298  

299  If the objective is found to contain only integer variables, set the 

300  fathoming discipline to strict. 

301  */ 

302  double maximumCost = 0.0 ; 

303  bool possibleMultiple = true ; 

304  int iColumn ; 

305  int numberColumns = getNumCols() ; 

306  for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) 

307  { if (upper[iColumn] > lower[iColumn]+1.0e8) 

308  { if (isInteger(iColumn)) 

309  maximumCost = CoinMax(maximumCost,fabs(objective[iColumn])) ; 

310  else if (objective[iColumn]) 

311  possibleMultiple = false ; } } 

312  setIntParam(CbcModel::CbcFathomDiscipline,possibleMultiple) ; 

313  /* 

314  If a nontrivial increment is possible, try and figure it out. We're looking 

315  for gcd(c<j>) for all c<j> that are coefficients of unfixed integer 

316  variables. Since the c<j> might not be integers, try and inflate them 

317  sufficiently that they look like integers (and we'll deflate the gcd 

318  later). 

319  

320  2520.0 is used as it is a nice multiple of 2,3,5,7 

321  */ 

322  if (possibleMultiple&&maximumCost) 

323  { int increment = 0 ; 

324  double multiplier = 2520.0 ; 

325  while (10.0*multiplier*maximumCost < 1.0e8) 

326  multiplier *= 10.0 ; 

327  

328  for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) 

329  { if (upper[iColumn] > lower[iColumn]+1.0e8) 

330  { if (isInteger(iColumn)&&objective[iColumn]) 

331  { double value = fabs(objective[iColumn])*multiplier ; 

332  int nearest = (int) floor(value+0.5) ; 

333  if (fabs(valuefloor(value+0.5)) > 1.0e8) 

334  { increment = 0 ; 

335  break ; } 

336  else if (!increment) 

337  { increment = nearest ; } 

338  else 

339  { increment = gcd(increment,nearest) ; } } } } 

340  /* 

341  If the increment beats the current value for objective change, install it. 

342  */ 

343  if (increment) 

344  { double value = increment ; 

345  double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ; 

346  value /= multiplier ; 

347  if (value*0.999 > cutoff) 

348  { messageHandler()>message(CBC_INTEGERINCREMENT, 

349  messages()) 

350  << value << CoinMessageEol ; 

351  setDblParam(CbcModel::CbcCutoffIncrement,value*0.999) ; } } } 

352  

353  return ; 

354  } 

355  

356  

357  /** 

358  \todo 

359  Normally, it looks like we enter here from command dispatch in the main 

360  routine, after calling the solver for an initial solution 

361  (CbcModel::initialSolve, which simply calls the solver's initialSolve 

362  routine.) The first thing we do is call resolve. Presumably there are 

363  circumstances where this is nontrivial? There's also a call from 

364  CbcModel::originalModel (tied up with integer presolve), which should be 

365  checked. 

366  

367  */ 

368  

369  /* 

370  The overall flow can be divided into three stages: 

371  * Prep: Check that the lp relaxation remains feasible at the root. If so, 

372  do all the setup for B&C. 

373  * Process the root node: Generate cuts, apply heuristics, and in general do 

374  the best we can to resolve the problem without B&C. 

375  * Do B&C search until we hit a limit or exhaust the search tree. 

376  

377  Keep in mind that in general there is no node in the search tree that 

378  corresponds to the active subproblem. The active subproblem is represented 

379  by the current state of the model, of the solver, and of the constraint 

380  system held by the solver. 

381  */ 

382  

383  void CbcModel::branchAndBound(int doStatistics) 

384  

385  { 

386  /* 

387  Capture a time stamp before we start. 

388  */ 

389  dblParam_[CbcStartSeconds] = CoinCpuTime(); 

390  strongInfo_[0]=0; 

391  strongInfo_[1]=0; 

392  strongInfo_[2]=0; 

393  numberStrongIterations_ = 0; 

394  // original solver (only set if preprocessing) 

395  OsiSolverInterface * originalSolver=NULL; 

396  int numberOriginalObjects=numberObjects_; 

397  CbcObject ** originalObject = NULL; 

398  // Set up strategies 

399  if (strategy_) { 

400  // May do preprocessing 

401  originalSolver = solver_; 

402  strategy_>setupOther(*this); 

403  if (strategy_>preProcessState()) { 

404  // preprocessing done 

405  if (strategy_>preProcessState()<0) { 

406  // infeasible 

407  handler_>message(CBC_INFEAS,messages_)<< CoinMessageEol ; 

408  status_ = 0 ; 

409  secondaryStatus_ = 1; 

410  originalContinuousObjective_ = COIN_DBL_MAX; 

411  return ; 

412  } else if (numberObjects_&&object_) { 

413  numberOriginalObjects=numberObjects_; 

414  // redo sequence 

415  numberIntegers_=0; 

416  int numberColumns = getNumCols(); 

417  int nOrig = originalSolver>getNumCols(); 

418  CglPreProcess * process = strategy_>process(); 

419  assert (process); 

420  const int * originalColumns = process>originalColumns(); 

421  // allow for cliques etc 

422  nOrig = CoinMax(nOrig,originalColumns[numberColumns1]+1); 

423  originalObject = object_; 

424  // object number or 1 

425  int * temp = new int[nOrig]; 

426  int iColumn; 

427  for (iColumn=0;iColumn<nOrig;iColumn++) 

428  temp[iColumn]=1; 

429  int iObject; 

430  int nNonInt=0; 

431  for (iObject=0;iObject<numberOriginalObjects;iObject++) { 

432  iColumn = originalObject[iObject]>columnNumber(); 

433  if (iColumn<0) { 

434  nNonInt++; 

435  } else { 

436  temp[iColumn]=iObject; 

437  } 

438  } 

439  int numberNewIntegers=0; 

440  int numberOldIntegers=0; 

441  int numberOldOther=0; 

442  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

443  int jColumn = originalColumns[iColumn]; 

444  if (temp[jColumn]>=0) { 

445  int iObject= temp[jColumn]; 

446  CbcSimpleInteger * obj = 

447  dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ; 

448  if (obj) 

449  numberOldIntegers++; 

450  else 

451  numberOldOther++; 

452  } else if (isInteger(iColumn)) { 

453  numberNewIntegers++; 

454  } 

455  } 

456  /* 

457  Allocate an array to hold the indices of the integer variables. 

458  Make a large enough array for all objects 

459  */ 

460  numberObjects_= numberNewIntegers+numberOldIntegers+numberOldOther+nNonInt; 

461  object_ = new CbcObject * [numberObjects_]; 

462  integerVariable_ = new int [numberNewIntegers+numberOldIntegers]; 

463  /* 

464  Walk the variables again, filling in the indices and creating objects for 

465  the integer variables. Initially, the objects hold the index and upper & 

466  lower bounds. 

467  */ 

468  numberIntegers_=0; 

469  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

470  int jColumn = originalColumns[iColumn]; 

471  if (temp[jColumn]>=0) { 

472  int iObject= temp[jColumn]; 

473  CbcSimpleInteger * obj = 

474  dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ; 

475  if (obj) { 

476  object_[numberIntegers_] = originalObject[iObject]>clone(); 

477  // redo ids etc 

478  object_[numberIntegers_]>redoSequenceEtc(this,numberColumns,originalColumns); 

479  integerVariable_[numberIntegers_++]=iColumn; 

480  } 

481  } else if (isInteger(iColumn)) { 

482  object_[numberIntegers_] = 

483  new CbcSimpleInteger(this,numberIntegers_,iColumn); 

484  integerVariable_[numberIntegers_++]=iColumn; 

485  } 

486  } 

487  numberObjects_=numberIntegers_; 

488  // Now append other column stuff 

489  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

490  int jColumn = originalColumns[iColumn]; 

491  if (temp[jColumn]>=0) { 

492  int iObject= temp[jColumn]; 

493  CbcSimpleInteger * obj = 

494  dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ; 

495  if (!obj) { 

496  object_[numberObjects_] = originalObject[iObject]>clone(); 

497  // redo ids etc 

498  object_[numberObjects_]>redoSequenceEtc(this,numberColumns,originalColumns); 

499  numberObjects_++; 

500  } 

501  } 

502  } 

503  // now append non column stuff 

504  for (iObject=0;iObject<numberOriginalObjects;iObject++) { 

505  iColumn = originalObject[iObject]>columnNumber(); 

506  if (iColumn<0) { 

507  object_[numberObjects_] = originalObject[iObject]>clone(); 

508  // redo ids etc 

509  object_[numberObjects_]>redoSequenceEtc(this,numberColumns,originalColumns); 

510  numberObjects_++; 

511  } 

512  } 

513  delete [] temp; 

514  if (!numberObjects_) 

515  handler_>message(CBC_NOINT,messages_) << CoinMessageEol ; 

516  } 

517  } else { 

518  //no preprocessing 

519  originalSolver=NULL; 

520  } 

521  strategy_>setupCutGenerators(*this); 

522  strategy_>setupHeuristics(*this); 

523  // Set strategy print level to models 

524  strategy_>setupPrinting(*this,handler_>logLevel()); 

525  } 

526  eventHappened_=false; 

527  CbcEventHandler *eventHandler = getEventHandler() ; 

528  

529  if (!nodeCompare_) 

530  nodeCompare_=new CbcCompareDefault();; 

531  // See if hot start wanted 

532  CbcCompareBase * saveCompare = NULL; 

533  if (hotstartSolution_) { 

534  if (strategy_&&strategy_>preProcessState()>0) { 

535  CglPreProcess * process = strategy_>process(); 

536  assert (process); 

537  int n = solver_>getNumCols(); 

538  const int * originalColumns = process>originalColumns(); 

539  // columns should be in order ... but 

540  double * tempS = new double[n]; 

541  for (int i=0;i<n;i++) { 

542  int iColumn = originalColumns[i]; 

543  tempS[i]=hotstartSolution_[iColumn]; 

544  } 

545  delete [] hotstartSolution_; 

546  hotstartSolution_=tempS; 

547  if (hotstartPriorities_) { 

548  int * tempP = new int [n]; 

549  for (int i=0;i<n;i++) { 

550  int iColumn = originalColumns[i]; 

551  tempP[i]=hotstartPriorities_[iColumn]; 

552  } 

553  delete [] hotstartPriorities_; 

554  hotstartPriorities_=tempP; 

555  } 

556  } 

557  saveCompare = nodeCompare_; 

558  // depth first 

559  nodeCompare_ = new CbcCompareDepth(); 

560  } 

561  if (!problemFeasibility_) 

562  problemFeasibility_=new CbcFeasibilityBase(); 

563  # ifdef CBC_DEBUG 

564  std::string problemName ; 

565  solver_>getStrParam(OsiProbName,problemName) ; 

566  printf("Problem name  %s\n",problemName.c_str()) ; 

567  solver_>setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ; 

568  # endif 

569  /* 

570  Assume we're done, and see if we're proven wrong. 

571  */ 

572  status_ = 0 ; 

573  secondaryStatus_ = 0; 

574  phase_=0; 

575  /* 

576  Scan the variables, noting the integer variables. Create an 

577  CbcSimpleInteger object for each integer variable. 

578  */ 

579  findIntegers(false) ; 

580  // If dynamic pseudo costs then do 

581  if (numberBeforeTrust_) 

582  convertToDynamic(); 

583  // Set up char array to say if integer 

584  delete [] integerInfo_; 

585  { 

586  int n = solver_>getNumCols(); 

587  integerInfo_ = new char [n]; 

588  for (int i=0;i<n;i++) { 

589  if (solver_>isInteger(i)) 

590  integerInfo_[i]=1; 

591  else 

592  integerInfo_[i]=0; 

593  } 

594  } 

595  

596  /* 

597  Ensure that objects on the lists of CbcObjects, heuristics, and cut 

598  generators attached to this model all refer to this model. 

599  */ 

600  synchronizeModel() ; 

601  assert (!solverCharacteristics_); 

602  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_>getAuxiliaryInfo()); 

603  if (solverCharacteristics) { 

604  solverCharacteristics_ = solverCharacteristics; 

605  } else { 

606  // replace in solver 

607  OsiBabSolver defaultC; 

608  solver_>setAuxiliaryInfo(&defaultC); 

609  solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_>getAuxiliaryInfo()); 

610  } 

611  solverCharacteristics_>setSolver(solver_); 

612  // Set so we can tell we are in initial phase in resolve 

613  continuousObjective_ = COIN_DBL_MAX ; 

614  /* 

615  Solve the relaxation. 

616  

617  Apparently there are circumstances where this will be nontrivial  i.e., 

618  we've done something since initialSolve that's trashed the solution to the 

619  continuous relaxation. 

620  */ 

621  bool feasible; 

622  // If NLP then we assume already solved outside branchAndbound 

623  if (!solverCharacteristics_>solverType()) { 

624  feasible=resolve(NULL,0) != 0 ; 

625  } else { 

626  // pick up given status 

627  feasible = (solver_>isProvenOptimal() && 

628  !solver_>isDualObjectiveLimitReached()) ; 

629  } 

630  if (problemFeasibility_>feasible(this,0)<0) { 

631  feasible=false; // pretend infeasible 

632  } 

633  /* 

634  If the linear relaxation of the root is infeasible, bail out now. Otherwise, 

635  continue with processing the root node. 

636  */ 

637  if (!feasible) 

638  { handler_>message(CBC_INFEAS,messages_)<< CoinMessageEol ; 

639  status_ = 0 ; 

640  secondaryStatus_ = 1; 

641  originalContinuousObjective_ = COIN_DBL_MAX; 

642  solverCharacteristics_ = NULL; 

643  return ; } 

644  // Save objective (just so user can access it) 

645  originalContinuousObjective_ = solver_>getObjValue(); 

646  bestPossibleObjective_=originalContinuousObjective_; 

647  sumChangeObjective1_=0.0; 

648  sumChangeObjective2_=0.0; 

649  /* 

650  OsiRowCutDebugger knows an optimal answer for a subset of MIP problems. 

651  Assuming it recognises the problem, when called upon it will check a cut to 

652  see if it cuts off the optimal answer. 

653  */ 

654  // If debugger exists set specialOptions_ bit 

655  if (solver_>getRowCutDebuggerAlways()) 

656  specialOptions_ = 1; 

657  

658  # ifdef CBC_DEBUG 

659  if ((specialOptions_&1)==0) 

660  solver_>activateRowCutDebugger(problemName.c_str()) ; 

661  if (solver_>getRowCutDebuggerAlways()) 

662  specialOptions_ = 1; 

663  # endif 

664  

665  /* 

666  Begin setup to process a feasible root node. 

667  */ 

668  bestObjective_ = CoinMin(bestObjective_,1.0e50) ; 

669  numberSolutions_ = 0 ; 

670  stateOfSearch_=0; 

671  numberHeuristicSolutions_ = 0 ; 

672  // Everything is minimization 

673  { 

674  // needed to sync cutoffs 

675  double value ; 

676  solver_>getDblParam(OsiDualObjectiveLimit,value) ; 

677  dblParam_[CbcCurrentCutoff]= value * solver_>getObjSense(); 

678  } 

679  double cutoff=getCutoff() ; 

680  double direction = solver_>getObjSense() ; 

681  dblParam_[CbcOptimizationDirection]=direction; 

682  if (cutoff < 1.0e20&&direction<0.0) 

683  messageHandler()>message(CBC_CUTOFF_WARNING1, 

684  messages()) 

685  << cutoff << cutoff << CoinMessageEol ; 

686  if (cutoff > bestObjective_) 

687  cutoff = bestObjective_ ; 

688  setCutoff(cutoff) ; 

689  /* 

690  We probably already have a current solution, but just in case ... 

691  */ 

692  int numberColumns = getNumCols() ; 

693  if (!currentSolution_) 

694  currentSolution_ = new double[numberColumns] ; 

695  testSolution_ = currentSolution_; 

696  /* Tell solver we are in Branch and Cut 

697  Could use last parameter for subtle differences */ 

698  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ; 

699  /* 

700  Create a copy of the solver, thus capturing the original (root node) 

701  constraint system (aka the continuous system). 

702  */ 

703  continuousSolver_ = solver_>clone() ; 

704  #ifdef CBC_USE_CLP 

705  OsiClpSolverInterface * clpSolver 

706  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

707  if (clpSolver) { 

708  ClpSimplex * clpSimplex = clpSolver>getModelPtr(); 

709  // take off names 

710  clpSimplex>dropNames(); 

711  } 

712  #endif 

713  

714  numberRowsAtContinuous_ = getNumRows() ; 

715  /* 

716  Check the objective to see if we can deduce a nontrivial increment. If 

717  it's better than the current value for CbcCutoffIncrement, it'll be 

718  installed. 

719  */ 

720  if(solverCharacteristics_>reducedCostsAccurate()) 

721  analyzeObjective() ; 

722  /* 

723  Set up for cut generation. addedCuts_ holds the cuts which are relevant for 

724  the active subproblem. whichGenerator will be used to record the generator 

725  that produced a given cut. 

726  */ 

727  maximumWhich_ = 1000 ; 

728  delete [] whichGenerator_; 

729  whichGenerator_ = new int[maximumWhich_] ; 

730  memset(whichGenerator_,0,maximumWhich_*sizeof(int)); 

731  int currentNumberCuts = 0 ; 

732  maximumNumberCuts_ = 0 ; 

733  currentNumberCuts_ = 0 ; 

734  delete [] addedCuts_ ; 

735  addedCuts_ = NULL ; 

736  /* 

737  Set up an empty heap and associated data structures to hold the live set 

738  (problems which require further exploration). 

739  */ 

740  tree_>setComparison(*nodeCompare_) ; 

741  /* 

742  Used to record the path from a node to the root of the search tree, so that 

743  we can then traverse from the root to the node when restoring a subproblem. 

744  */ 

745  maximumDepth_ = 10 ; 

746  delete [] walkback_ ; 

747  walkback_ = new CbcNodeInfo * [maximumDepth_] ; 

748  /* 

749  Used to generate bound edits for CbcPartialNodeInfo. 

750  */ 

751  double * lowerBefore = new double [numberColumns] ; 

752  double * upperBefore = new double [numberColumns] ; 

753  /* 

754  

755  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy 

756  lifting. It will iterate a generate/reoptimise loop (including reduced cost 

757  fixing) until no cuts are generated, the change in objective falls off, or 

758  the limit on the number of rounds of cut generation is exceeded. 

759  

760  At the end of all this, any cuts will be recorded in cuts and also 

761  installed in the solver's constraint system. We'll have reoptimised, and 

762  removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been 

763  adjusted accordingly). 

764  

765  Tell cut generators they can be a bit more aggressive at root node 

766  

767  TODO: Why don't we make a copy of the solution after solveWithCuts? 

768  TODO: If numberUnsatisfied == 0, don't we have a solution? 

769  */ 

770  phase_=1; 

771  int iCutGenerator; 

772  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) { 

773  CglCutGenerator * generator = generator_[iCutGenerator]>generator(); 

774  generator>setAggressiveness(generator>getAggressiveness()+100); 

775  } 

776  OsiCuts cuts ; 

777  int anyAction = 1 ; 

778  numberOldActiveCuts_ = 0 ; 

779  numberNewCuts_ = 0 ; 

780  // Array to mark solution 

781  delete [] usedInSolution_; 

782  usedInSolution_ = new int[numberColumns]; 

783  CoinZeroN(usedInSolution_,numberColumns); 

784  /* 

785  For printing totals and for CbcNode (numberNodes_) 

786  */ 

787  numberIterations_ = 0 ; 

788  numberNodes_ = 0 ; 

789  numberNodes2_ = 0 ; 

790  maximumStatistics_=0; 

791  statistics_ = NULL; 

792  // Do on switch 

793  if (doStatistics) { 

794  maximumStatistics_=10000; 

795  statistics_ = new CbcStatistics * [maximumStatistics_]; 

796  memset(statistics_,0,maximumStatistics_*sizeof(CbcStatistics *)); 

797  } 

798  

799  { int iObject ; 

800  int preferredWay ; 

801  int numberUnsatisfied = 0 ; 

802  memcpy(currentSolution_,solver_>getColSolution(), 

803  numberColumns*sizeof(double)) ; 

804  

805  for (iObject = 0 ; iObject < numberObjects_ ; iObject++) 

806  { double infeasibility = 

807  object_[iObject]>infeasibility(preferredWay) ; 

808  if (infeasibility ) numberUnsatisfied++ ; } 

809  // replace solverType 

810  if(solverCharacteristics_>tryCuts()) { 

811  if (numberUnsatisfied) { 

812  feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_, 

813  NULL); 

814  } else if (solverCharacteristics_>solutionAddsCuts()) { 

815  // may generate cuts and turn the solution 

816  //to an infeasible one 

817  feasible = solveWithCuts(cuts, 1, 

818  NULL); 

819  #if 0 

820  currentNumberCuts_ = cuts.sizeRowCuts(); 

821  if (currentNumberCuts_ >= maximumNumberCuts_) { 

822  maximumNumberCuts_ = currentNumberCuts; 

823  delete [] addedCuts_; 

824  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

825  } 

826  #endif 

827  } 

828  } 

829  // check extra info on feasibility 

830  if (!solverCharacteristics_>mipFeasible()) 

831  feasible = false; 

832  } 

833  // make cut generators less aggressive 

834  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) { 

835  CglCutGenerator * generator = generator_[iCutGenerator]>generator(); 

836  generator>setAggressiveness(generator>getAggressiveness()100); 

837  } 

838  currentNumberCuts_ = numberNewCuts_ ; 

839  /* 

840  We've taken the continuous relaxation as far as we can. Time to branch. 

841  The first order of business is to actually create a node. chooseBranch 

842  currently uses strong branching to evaluate branch object candidates, 

843  unless forced back to simple branching. If chooseBranch concludes that a 

844  branching candidate is monotone (anyAction == 1) or infeasible (anyAction 

845  == 2) when forced to integer values, it returns here immediately. 

846  

847  Monotone variables trigger a call to resolve(). If the problem remains 

848  feasible, try again to choose a branching variable. At the end of the loop, 

849  resolved == true indicates that some variables were fixed. 

850  

851  Loss of feasibility will result in the deletion of newNode. 

852  */ 

853  

854  bool resolved = false ; 

855  CbcNode *newNode = NULL ; 

856  numberFixedAtRoot_=0; 

857  numberFixedNow_=0; 

858  int numberIterationsAtContinuous = numberIterations_; 

859  //solverCharacteristics_>setSolver(solver_); 

860  if (feasible) { 

861  newNode = new CbcNode ; 

862  // Set objective value (not so obvious if NLP etc) 

863  setObjectiveValue(newNode,NULL); 

864  anyAction = 1 ; 

865  // To make depth available we may need a fake node 

866  CbcNode fakeNode; 

867  if (!currentNode_) { 

868  // Not true if sub trees assert (!numberNodes_); 

869  currentNode_=&fakeNode; 

870  } 

871  phase_=3; 

872  // only allow 1000 passes 

873  int numberPassesLeft=1000; 

874  // This is first crude step 

875  if (numberAnalyzeIterations_) { 

876  delete [] analyzeResults_; 

877  analyzeResults_ = new double [4*numberIntegers_]; 

878  numberFixedAtRoot_=newNode>analyze(this,analyzeResults_); 

879  if (numberFixedAtRoot_>0) { 

880  printf("%d fixed by analysis\n",numberFixedAtRoot_); 

881  setPointers(solver_); 

882  numberFixedNow_ = numberFixedAtRoot_; 

883  } else if (numberFixedAtRoot_<0) { 

884  printf("analysis found to be infeasible\n"); 

885  anyAction=2; 

886  delete newNode ; 

887  newNode = NULL ; 

888  feasible = false ; 

889  } 

890  } 

891  while (anyAction == 1) 

892  { 

893  // Set objective value (not so obvious if NLP etc) 

894  setObjectiveValue(newNode,NULL); 

895  if (numberBeforeTrust_==0 ) { 

896  anyAction = newNode>chooseBranch(this,NULL,numberPassesLeft) ; 

897  } else { 

898  OsiSolverBranch * branches=NULL; 

899  anyAction = newNode>chooseDynamicBranch(this,NULL,branches,numberPassesLeft) ; 

900  if (anyAction==3) 

901  anyAction = newNode>chooseBranch(this,NULL,numberPassesLeft) ; // dynamic did nothing 

902  } 

903  numberPassesLeft; 

904  if (anyAction == 1) 

905  { feasible = resolve(NULL,10) != 0 ; 

906  if (problemFeasibility_>feasible(this,0)<0) { 

907  feasible=false; // pretend infeasible 

908  } 

909  reducedCostFix() ; 

910  resolved = true ; 

911  # ifdef CBC_DEBUG 

912  printf("Resolve (root) as something fixed, Obj value %g %d rows\n", 

913  solver_>getObjValue(), 

914  solver_>getNumRows()) ; 

915  # endif 

916  if (!feasible) anyAction = 2 ; } 

917  if (anyAction == 2newNode>objectiveValue() >= cutoff) 

918  { delete newNode ; 

919  newNode = NULL ; 

920  feasible = false ; } } } 

921  /* 

922  At this point, the root subproblem is infeasible or fathomed by bound 

923  (newNode == NULL), or we're live with an objective value that satisfies the 

924  current objective cutoff. 

925  */ 

926  assert (!newNode  newNode>objectiveValue() <= cutoff) ; 

927  // Save address of root node as we don't want to delete it 

928  CbcNode * rootNode = newNode; 

929  /* 

930  The common case is that the lp relaxation is feasible but doesn't satisfy 

931  integrality (i.e., newNode>variable() >= 0, indicating we've been able to 

932  select a branching variable). Remove any cuts that have gone slack due to 

933  forcing monotone variables. Then tack on an CbcFullNodeInfo object and full 

934  basis (via createInfo()) and stash the new cuts in the nodeInfo (via 

935  addCuts()). If, by some miracle, we have an integral solution at the root 

936  (newNode>variable() < 0), takeOffCuts() will ensure that the solver holds 

937  a valid solution for use by setBestSolution(). 

938  */ 

939  CoinWarmStartBasis *lastws = 0 ; 

940  if (feasible && newNode>variable() >= 0) 

941  { if (resolved) 

942  { bool needValidSolution = (newNode>variable() < 0) ; 

943  takeOffCuts(cuts,needValidSolution,NULL) ; 

944  # ifdef CHECK_CUT_COUNTS 

945  { printf("Number of rows after chooseBranch fix (root)" 

946  "(active only) %d\n", 

947  numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ; 

948  const CoinWarmStartBasis* debugws = 

949  dynamic_cast <const CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

950  debugws>print() ; 

951  delete debugws ; } 

952  # endif 

953  } 

954  newNode>createInfo(this,NULL,NULL,NULL,NULL,0,0) ; 

955  newNode>nodeInfo()>addCuts(cuts, 

956  newNode>numberBranches(),whichGenerator_) ; 

957  if (lastws) delete lastws ; 

958  lastws = dynamic_cast<CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

959  } 

960  /* 

961  Continuous data to be used later 

962  */ 

963  continuousObjective_ = solver_>getObjValue()*solver_>getObjSense(); 

964  continuousInfeasibilities_ = 0 ; 

965  if (newNode) 

966  { continuousObjective_ = newNode>objectiveValue() ; 

967  delete [] continuousSolution_; 

968  continuousSolution_ = CoinCopyOfArray(solver_>getColSolution(), 

969  numberColumns); 

970  continuousInfeasibilities_ = newNode>numberUnsatisfied() ; } 

971  /* 

972  Bound may have changed so reset in objects 

973  */ 

974  { int i ; 

975  for (i = 0;i < numberObjects_;i++) 

976  object_[i]>resetBounds() ; } 

977  stoppedOnGap_ = false ; 

978  /* 

979  Feasible? Then we should have either a live node prepped for future 

980  expansion (indicated by variable() >= 0), or (miracle of miracles) an 

981  integral solution at the root node. 

982  

983  initializeInfo sets the reference counts in the nodeInfo object. Since 

984  this node is still live, push it onto the heap that holds the live set. 

985  */ 

986  double bestValue = 0.0 ; 

987  if (newNode) { 

988  bestValue = newNode>objectiveValue(); 

989  if (newNode>variable() >= 0) { 

990  newNode>initializeInfo() ; 

991  tree_>push(newNode) ; 

992  if (statistics_) { 

993  if (numberNodes2_==maximumStatistics_) { 

994  maximumStatistics_ = 2*maximumStatistics_; 

995  CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_]; 

996  memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *)); 

997  memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *)); 

998  delete [] statistics_; 

999  statistics_=temp; 

1000  } 

1001  assert (!statistics_[numberNodes2_]); 

1002  statistics_[numberNodes2_]=new CbcStatistics(newNode); 

1003  } 

1004  numberNodes2_++; 

1005  # ifdef CHECK_NODE 

1006  printf("Node %x on tree\n",newNode) ; 

1007  # endif 

1008  } else { 

1009  // continuous is integer 

1010  double objectiveValue = newNode>objectiveValue(); 

1011  setBestSolution(CBC_SOLUTION,objectiveValue, 

1012  solver_>getColSolution()) ; 

1013  delete newNode ; 

1014  newNode = NULL ; 

1015  } 

1016  } 

1017  

1018  if (printFrequency_ <= 0) { 

1019  printFrequency_ = 1000 ; 

1020  if (getNumCols() > 2000) 

1021  printFrequency_ = 100 ; 

1022  } 

1023  /* 

1024  It is possible that strong branching fixes one variable and then the code goes round 

1025  again and again. This can take too long. So we need to warn user  just once. 

1026  */ 

1027  numberLongStrong_=0; 

1028  /* 

1029  At last, the actual branchandcut search loop, which will iterate until 

1030  the live set is empty or we hit some limit (integrality gap, time, node 

1031  count, etc.). The overall flow is to rebuild a subproblem, reoptimise using 

1032  solveWithCuts(), choose a branching pattern with chooseBranch(), and finally 

1033  add the node to the live set. 

1034  

1035  The first action is to winnow the live set to remove nodes which are worse 

1036  than the current objective cutoff. 

1037  */ 

1038  while (!tree_>empty()) { 

1039  #ifdef BONMIN 

1040  assert(!solverCharacteristics_>solutionAddsCuts()  solverCharacteristics_>mipFeasible()); 

1041  #endif 

1042  if (cutoff > getCutoff()) { 

1043  double newCutoff = getCutoff(); 

1044  if (analyzeResults_) { 

1045  // see if we could fix any (more) 

1046  int n=0; 

1047  double * newLower = analyzeResults_; 

1048  double * objLower = newLower+numberIntegers_; 

1049  double * newUpper = objLower+numberIntegers_; 

1050  double * objUpper = newUpper+numberIntegers_; 

1051  for (int i=0;i<numberIntegers_;i++) { 

1052  if (objLower[i]>newCutoff) { 

1053  n++; 

1054  if (objUpper[i]>newCutoff) { 

1055  newCutoff = COIN_DBL_MAX; 

1056  break; 

1057  } 

1058  } else if (objUpper[i]>newCutoff) { 

1059  n++; 

1060  } 

1061  } 

1062  if (newCutoff==COIN_DBL_MAX) { 

1063  printf("Root analysis says finished\n"); 

1064  } else if (n>numberFixedNow_) { 

1065  printf("%d more fixed by analysis  now %d\n",nnumberFixedNow_,n); 

1066  numberFixedNow_=n; 

1067  } 

1068  } 

1069  if (eventHandler) { 

1070  if (!eventHandler>event(CbcEventHandler::solution)) { 

1071  eventHappened_=true; // exit 

1072  } 

1073  } 

1074  // Do from deepest 

1075  tree_>cleanTree(this, newCutoff,bestPossibleObjective_) ; 

1076  nodeCompare_>newSolution(this) ; 

1077  nodeCompare_>newSolution(this,continuousObjective_, 

1078  continuousInfeasibilities_) ; 

1079  tree_>setComparison(*nodeCompare_) ; 

1080  if (tree_>empty()) 

1081  break; // finished 

1082  } 

1083  cutoff = getCutoff() ; 

1084  /* 

1085  Periodic activities: Opportunities to 

1086  + tweak the nodeCompare criteria, 

1087  + check if we've closed the integrality gap enough to quit, 

1088  + print a summary line to let the user know we're working 

1089  */ 

1090  if ((numberNodes_%1000) == 0) { 

1091  bool redoTree=nodeCompare_>every1000Nodes(this, numberNodes_) ; 

1092  // redo tree if wanted 

1093  if (redoTree) 

1094  tree_>setComparison(*nodeCompare_) ; 

1095  } 

1096  if (saveCompare&&!hotstartSolution_) { 

1097  // hotstart switched off 

1098  delete nodeCompare_; // off depth first 

1099  nodeCompare_=saveCompare; 

1100  saveCompare=NULL; 

1101  // redo tree 

1102  tree_>setComparison(*nodeCompare_) ; 

1103  } 

1104  if ((numberNodes_%printFrequency_) == 0) { 

1105  int j ; 

1106  int nNodes = tree_>size() ; 

1107  bestPossibleObjective_ = 1.0e100 ; 

1108  for (j = 0;j < nNodes;j++) { 

1109  CbcNode * node = tree_>nodePointer(j) ; 

1110  if (node&&node>objectiveValue() < bestPossibleObjective_) 

1111  bestPossibleObjective_ = node>objectiveValue() ; 

1112  } 

1113  messageHandler()>message(CBC_STATUS,messages()) 

1114  << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_ 

1115  <<getCurrentSeconds() 

1116  << CoinMessageEol ; 

1117  if (!eventHandler>event(CbcEventHandler::treeStatus)) { 

1118  eventHappened_=true; // exit 

1119  } 

1120  } 

1121  // If no solution but many nodes  signal change in strategy 

1122  if (numberNodes_>2*numberObjects_+1000&&stateOfSearch_!=2) 

1123  stateOfSearch_=3; 

1124  // See if can stop on gap 

1125  double testGap = CoinMax(dblParam_[CbcAllowableGap], 

1126  CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_)) 

1127  *dblParam_[CbcAllowableFractionGap]); 

1128  if (bestObjective_bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) { 

1129  stoppedOnGap_ = true ; 

1130  } 

1131  

1132  # ifdef CHECK_NODE_FULL 

1133  verifyTreeNodes(tree_,*this) ; 

1134  # endif 

1135  # ifdef CHECK_CUT_COUNTS 

1136  verifyCutCounts(tree_,*this) ; 

1137  # endif 

1138  

1139  /* 

1140  Now we come to the meat of the loop. To create the active subproblem, we'll 

1141  pop the most promising node in the live set, rebuild the subproblem it 

1142  represents, and then execute the current arm of the branch to create the 

1143  active subproblem. 

1144  */ 

1145  CbcNode *node = tree_>bestNode(cutoff) ; 

1146  // Possible one on tree worse than cutoff 

1147  if (!node) 

1148  continue; 

1149  currentNode_=node; // so can be accessed elsewhere 

1150  #ifdef CBC_DEBUG 

1151  printf("%d unsat, way %d, obj %g est %g\n", 

1152  node>numberUnsatisfied(),node>way(),node>objectiveValue(), 

1153  node>guessedObjectiveValue()); 

1154  #endif 

1155  // Save clone in branching decision 

1156  if(branchingMethod_) 

1157  branchingMethod_>saveBranchingObject(node>modifiableBranchingObject()); 

1158  bool nodeOnTree=false; // Node has been popped 

1159  // Say not on optimal path 

1160  bool onOptimalPath=false; 

1161  # ifdef CHECK_NODE 

1162  /* 

1163  WARNING: The use of integerVariable_[*] here will break as soon as the 

1164  branching object is something other than an integer variable. 

1165  This needs some thought. 

1166  */ 

1167  printf("Node %x popped from tree  %d left, %d count\n",node, 

1168  node>nodeInfo()>numberBranchesLeft(), 

1169  node>nodeInfo()>numberPointingToThis()) ; 

1170  printf("\tdepth = %d, z = %g, unsat = %d, var = %d.\n", 

1171  node>depth(),node>objectiveValue(), 

1172  node>numberUnsatisfied(), 

1173  integerVariable_[node>variable()]) ; 

1174  # endif 

1175  

1176  /* 

1177  Rebuild the subproblem for this node: Call addCuts() to adjust the model 

1178  to recreate the subproblem for this node (set proper variable bounds, add 

1179  cuts, create a basis). This may result in the problem being fathomed by 

1180  bound or infeasibility. Returns 1 if node is fathomed. 

1181  Execute the current arm of the branch: If the problem survives, save the 

1182  resulting variable bounds and call branch() to modify variable bounds 

1183  according to the current arm of the branching object. If we're processing 

1184  the final arm of the branching object, flag the node for removal from the 

1185  live set. 

1186  */ 

1187  CbcNodeInfo * nodeInfo = node>nodeInfo() ; 

1188  newNode = NULL ; 

1189  if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_)) 

1190  { int i ; 

1191  const double * lower = getColLower() ; 

1192  const double * upper = getColUpper() ; 

1193  for (i = 0 ; i < numberColumns ; i++) 

1194  { lowerBefore[i]= lower[i] ; 

1195  upperBefore[i]= upper[i] ; } 

1196  bool deleteNode ; 

1197  if (node>branch()) 

1198  { 

1199  // set nodenumber correctly 

1200  node>nodeInfo()>setNodeNumber(numberNodes2_); 

1201  tree_>push(node) ; 

1202  if (statistics_) { 

1203  if (numberNodes2_==maximumStatistics_) { 

1204  maximumStatistics_ = 2*maximumStatistics_; 

1205  CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_]; 

1206  memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *)); 

1207  memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *)); 

1208  delete [] statistics_; 

1209  statistics_=temp; 

1210  } 

1211  assert (!statistics_[numberNodes2_]); 

1212  statistics_[numberNodes2_]=new CbcStatistics(node); 

1213  } 

1214  numberNodes2_++; 

1215  nodeOnTree=true; // back on tree 

1216  deleteNode = false ; 

1217  # ifdef CHECK_NODE 

1218  printf("Node %x pushed back on tree  %d left, %d count\n",node, 

1219  nodeInfo>numberBranchesLeft(), 

1220  nodeInfo>numberPointingToThis()) ; 

1221  # endif 

1222  } 

1223  else 

1224  { deleteNode = true ; 

1225  if (!nodeInfo>numberBranchesLeft()) 

1226  nodeInfo>allBranchesGone(); // can clean up 

1227  } 

1228  

1229  if ((specialOptions_&1)!=0) { 

1230  /* 

1231  This doesn't work as intended  getRowCutDebugger will return null 

1232  unless the current feasible solution region includes the optimal solution 

1233  that RowCutDebugger knows. There's no way to tell inactive from off the 

1234  optimal path. 

1235  */ 

1236  const OsiRowCutDebugger *debugger = solver_>getRowCutDebugger() ; 

1237  if (debugger) { 

1238  onOptimalPath=true; 

1239  printf("On optimal path\n") ; 

1240  } 

1241  } 

1242  /* 

1243  Reoptimize, possibly generating cuts and/or using heuristics to find 

1244  solutions. Cut reference counts are unaffected unless we lose feasibility, 

1245  in which case solveWithCuts() will make the adjustment. 

1246  */ 

1247  phase_=2; 

1248  cuts = OsiCuts() ; 

1249  currentNumberCuts = solver_>getNumRows()numberRowsAtContinuous_ ; 

1250  int saveNumber = numberIterations_; 

1251  if(solverCharacteristics_>solutionAddsCuts()) { 

1252  int returnCode=resolve(node ? node>nodeInfo() : NULL,1); 

1253  feasible = returnCode != 0; 

1254  if (feasible) { 

1255  int iObject ; 

1256  int preferredWay ; 

1257  int numberUnsatisfied = 0 ; 

1258  memcpy(currentSolution_,solver_>getColSolution(), 

1259  numberColumns*sizeof(double)) ; 

1260  

1261  for (iObject = 0 ; iObject < numberObjects_ ; iObject++) { 

1262  double infeasibility = 

1263  object_[iObject]>infeasibility(preferredWay) ; 

1264  if (infeasibility ) numberUnsatisfied++ ; 

1265  } 

1266  if (returnCode>0) { 

1267  if (numberUnsatisfied) { 

1268  feasible = solveWithCuts(cuts,maximumCutPasses_,node); 

1269  } else { 

1270  // may generate cuts and turn the solution 

1271  //to an infeasible one 

1272  feasible = solveWithCuts(cuts, 1, 

1273  node); 

1274  #if 0 

1275  currentNumberCuts_ = cuts.sizeRowCuts(); 

1276  if (currentNumberCuts_ >= maximumNumberCuts_) { 

1277  maximumNumberCuts_ = currentNumberCuts; 

1278  delete [] addedCuts_; 

1279  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

1280  } 

1281  #endif 

1282  } 

1283  } 

1284  // check extra info on feasibility 

1285  if (!solverCharacteristics_>mipFeasible()) 

1286  feasible = false; 

1287  } 

1288  } else { 

1289  // normal 

1290  feasible = solveWithCuts(cuts,maximumCutPasses_,node); 

1291  } 

1292  if ((specialOptions_&1)!=0&&onOptimalPath) { 

1293  const OsiRowCutDebugger *debugger = solver_>getRowCutDebugger() ; 

1294  assert (debugger) ; 

1295  } 

1296  if (statistics_) { 

1297  assert (numberNodes2_); 

1298  assert (statistics_[numberNodes2_1]); 

1299  assert (statistics_[numberNodes2_1]>node()==numberNodes2_1); 

1300  statistics_[numberNodes2_1]>endOfBranch(numberIterations_saveNumber, 

1301  feasible ? solver_>getObjValue() 

1302  : COIN_DBL_MAX); 

1303  } 

1304  /* 

1305  Check for abort on limits: node count, solution count, time, integrality gap. 

1306  */ 

1307  double totalTime = CoinCpuTime()dblParam_[CbcStartSeconds] ; 

1308  if (numberNodes_ < intParam_[CbcMaxNumNode] && 

1309  numberSolutions_ < intParam_[CbcMaxNumSol] && 

1310  totalTime < dblParam_[CbcMaximumSeconds] && 

1311  !stoppedOnGap_&&!eventHappened_) 

1312  { 

1313  /* 

1314  Are we still feasible? If so, create a node and do the work to attach a 

1315  branching object, reoptimising as needed if chooseBranch() identifies 

1316  monotone objects. 

1317  

1318  Finally, attach a partial nodeInfo object and store away any cuts that we 

1319  created back in solveWithCuts. addCuts() will initialise the reference 

1320  counts for these new cuts. 

1321  

1322  TODO: (lh) I'm confused. We create a nodeInfo without checking whether we 

1323  have a solution or not. Then we use numberUnsatisfied() to decide 

1324  whether to stash the cuts and bump reference counts. Other places we 

1325  use variable() (i.e., presence of a branching variable). Equivalent? 

1326  */ 

1327  if (onOptimalPath) { 

1328  if (!feasible) { 

1329  printf("infeas2\n"); 

1330  solver_>writeMps("infeas"); 

1331  CoinWarmStartBasis *slack = 

1332  dynamic_cast<CoinWarmStartBasis *>(solver_>getEmptyWarmStart()) ; 

1333  solver_>setWarmStart(slack); 

1334  delete slack ; 

1335  solver_>setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ; 

1336  solver_>initialSolve(); 

1337  assert (!solver_>isProvenOptimal()); 

1338  } 

1339  assert (feasible); 

1340  } 

1341  bool checkingNode=false; 

1342  if (feasible) { 

1343  newNode = new CbcNode ; 

1344  // Set objective value (not so obvious if NLP etc) 

1345  setObjectiveValue(newNode,node); 

1346  anyAction =1 ; 

1347  resolved = false ; 

1348  if (newNode>objectiveValue() >= getCutoff()) 

1349  anyAction=2; 

1350  // only allow twenty passes 

1351  int numberPassesLeft=20; 

1352  checkingNode=true; 

1353  while (anyAction == 1) 

1354  { 

1355  // Set objective value (not so obvious if NLP etc) 

1356  setObjectiveValue(newNode,node); 

1357  if (numberBeforeTrust_==0 ) { 

1358  anyAction = newNode>chooseBranch(this,node,numberPassesLeft) ; 

1359  } else { 

1360  OsiSolverBranch * branches=NULL; 

1361  anyAction = newNode>chooseDynamicBranch(this,node,branches,numberPassesLeft) ; 

1362  if (anyAction==3) 

1363  anyAction = newNode>chooseBranch(this,node,numberPassesLeft) ; // dynamic did nothing 

1364  } 

1365  if (solverCharacteristics_ && 

1366  solverCharacteristics_>solutionAddsCuts() && // we are in some OA based bab 

1367  feasible && (newNode>numberUnsatisfied()==0) //solution has become integer feasible during strong branching 

1368  ) 

1369  { 

1370  //in the present case we need to check here integer infeasibility if the node is not fathomed we will have to do the loop 

1371  // again 

1372  std::cout<<solver_<<std::endl; 

1373  solver_>resolve(); 

1374  double objval = solver_>getObjValue(); 

1375  setBestSolution(CBC_SOLUTION, objval, 

1376  solver_>getColSolution()) ; 

1377  lastHeuristic_ = NULL; 

1378  int easy=2; 

1379  if (!solverCharacteristics_>mipFeasible())//did we prove that the node could be pruned? 

1380  feasible = false; 

1381  // Reset the bound now 

1382  solverCharacteristics_>setMipBound(COIN_DBL_MAX); 

1383  

1384  

1385  //numberPassesLeft++; 

1386  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ; 

1387  feasible &= resolve(node ? node>nodeInfo() : NULL,11) != 0 ; 

1388  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ; 

1389  resolved = true ; 

1390  if (problemFeasibility_>feasible(this,0)<0) { 

1391  feasible=false; // pretend infeasible 

1392  } 

1393  if(feasible) 

1394  anyAction = 1; 

1395  else 

1396  anyAction = 2; 

1397  } 

1398  /* 

1399  Yep, false positives for sure. And no easy way to distinguish honest 

1400  infeasibility from `found a solution and tightened objective target.' 

1401  

1402  if (onOptimalPath) 

1403  assert (anyAction!=2); // can be useful but gives false positives on strong 

1404  */ 

1405  numberPassesLeft; 

1406  if (numberPassesLeft<=1) { 

1407  if (!numberLongStrong_) 

1408  messageHandler()>message(CBC_WARNING_STRONG, 

1409  messages()) << CoinMessageEol ; 

1410  numberLongStrong_++; 

1411  } 

1412  if (anyAction == 1) 

1413  { 

1414  // can do quick optimality check 

1415  int easy=2; 

1416  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ; 

1417  feasible = resolve(node ? node>nodeInfo() : NULL,11) != 0 ; 

1418  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ; 

1419  resolved = true ; 

1420  if (problemFeasibility_>feasible(this,0)<0) { 

1421  feasible=false; // pretend infeasible 

1422  } 

1423  if (feasible) 

1424  { 

1425  // Set objective value (not so obvious if NLP etc) 

1426  setObjectiveValue(newNode,node); 

1427  reducedCostFix() ; 

1428  if (newNode>objectiveValue() >= getCutoff()) 

1429  anyAction=2; 

1430  } 

1431  else 

1432  { anyAction = 2 ; } } } 

1433  if (anyAction >= 0) 

1434  { if (resolved) 

1435  { bool needValidSolution = (newNode>variable() < 0) ; 

1436  takeOffCuts(cuts,needValidSolution,NULL) ; 

1437  # ifdef CHECK_CUT_COUNTS 

1438  { printf("Number of rows after chooseBranch fix (node)" 

1439  "(active only) %d\n", 

1440  numberRowsAtContinuous_+numberNewCuts_+ 

1441  numberOldActiveCuts_) ; 

1442  const CoinWarmStartBasis* debugws = 

1443  dynamic_cast<const CoinWarmStartBasis*> 

1444  (solver_>getWarmStart()) ; 

1445  debugws>print() ; 

1446  delete debugws ; } 

1447  # endif 

1448  } 

1449  newNode>createInfo(this,node,lastws,lowerBefore,upperBefore, 

1450  numberOldActiveCuts_,numberNewCuts_) ; 

1451  if (newNode>numberUnsatisfied()) { 

1452  newNode>initializeInfo() ; 

1453  newNode>nodeInfo()>addCuts(cuts,newNode>numberBranches(), 

1454  whichGenerator_) ; } } } 

1455  else { 

1456  anyAction = 2 ; 

1457  // Reset bound anyway (no harm if not odd) 

1458  solverCharacteristics_>setMipBound(COIN_DBL_MAX); 

1459  } 

1460  // May have slipped through i.e. anyAction == 0 and objective above cutoff 

1461  // I think this will screw up cut reference counts if executed. 

1462  // We executed addCuts just above. (lh) 

1463  if ( anyAction >=0 ) { 

1464  assert (newNode); 

1465  if (newNode>objectiveValue() >= getCutoff()) 

1466  anyAction = 2; // say bad after all 

1467  } 

1468  /* 

1469  If we end up infeasible, we can delete the new node immediately. Since this 

1470  node won't be needing the cuts we collected, decrement the reference counts. 

1471  If we are feasible, then we'll be placing this node into the live set, so 

1472  increment the reference count in the current (parent) nodeInfo. 

1473  */ 

1474  if (anyAction == 2) 

1475  { delete newNode ; 

1476  newNode = NULL ; 

1477  // say strong doing well 

1478  if (checkingNode) 

1479  setSpecialOptions(specialOptions_8); 

1480  for (i = 0 ; i < currentNumberCuts_ ; i++) 

1481  { if (addedCuts_[i]) 

1482  { if (!addedCuts_[i]>decrement(1)) 

1483  delete addedCuts_[i] ; } } } 

1484  else 

1485  { nodeInfo>increment() ; 

1486  if ((numberNodes_%20)==0) { 

1487  // say strong not doing as well 

1488  setSpecialOptions(specialOptions_&~8); 

1489  } 

1490  } 

1491  /* 

1492  At this point, there are three possibilities: 

1493  * newNode is live and will require further branching to resolve 

1494  (variable() >= 0). Increment the cut reference counts by 

1495  numberBranches() to allow for use by children of this node, and 

1496  decrement by 1 because we've executed one arm of the branch of our 

1497  parent (consuming one reference). Before we push newNode onto the 

1498  search tree, try for a heuristic solution. 

1499  * We have a solution, in which case newNode is nonnull but we have no 

1500  branching variable. Decrement the cut counts and save the solution. 

1501  * The node was found to be infeasible, in which case it's already been 

1502  deleted, and newNode is null. 

1503  */ 

1504  if (!eventHandler>event(CbcEventHandler::node)) { 

1505  eventHappened_=true; // exit 

1506  } 

1507  assert (!newNode  newNode>objectiveValue() <= getCutoff()) ; 

1508  if (statistics_) { 

1509  assert (numberNodes2_); 

1510  assert (statistics_[numberNodes2_1]); 

1511  assert (statistics_[numberNodes2_1]>node()==numberNodes2_1); 

1512  if (newNode) 

1513  statistics_[numberNodes2_1]>updateInfeasibility(newNode>numberUnsatisfied()); 

1514  else 

1515  statistics_[numberNodes2_1]>sayInfeasible(); 

1516  } 

1517  if (newNode) 

1518  { if (newNode>variable() >= 0) 

1519  { handler_>message(CBC_BRANCH,messages_) 

1520  << numberNodes_<< newNode>objectiveValue() 

1521  << newNode>numberUnsatisfied()<< newNode>depth() 

1522  << CoinMessageEol ; 

1523  // Increment cut counts (taking off current) 

1524  int numberLeft = newNode>numberBranches() ; 

1525  for (i = 0;i < currentNumberCuts_;i++) 

1526  { if (addedCuts_[i]) 

1527  { 

1528  # ifdef CHECK_CUT_COUNTS 

1529  printf("Count on cut %x increased by %d\n",addedCuts_[i], 

1530  numberLeft1) ; 

1531  # endif 

1532  addedCuts_[i]>increment(numberLeft1) ; } } 

1533  

1534  double estValue = newNode>guessedObjectiveValue() ; 

1535  int found = 1 ; 

1536  // no  overhead on small problems solver_>resolve() ; // double check current optimal 

1537  // assert (!solver_>getIterationCount()); 

1538  double * newSolution = new double [numberColumns] ; 

1539  double heurValue = getCutoff() ; 

1540  int iHeur ; 

1541  for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++) 

1542  { double saveValue = heurValue ; 

1543  int ifSol = heuristic_[iHeur]>solution(heurValue,newSolution) ; 

1544  if (ifSol > 0) { 

1545  // new solution found 

1546  found = iHeur ; 

1547  incrementUsed(newSolution); 

1548  } 

1549  else 

1550  if (ifSol < 0) // just returning an estimate 

1551  { estValue = CoinMin(heurValue,estValue) ; 

1552  heurValue = saveValue ; } } 

1553  if (found >= 0) { 

1554  setBestSolution(CBC_ROUNDING,heurValue,newSolution) ; 

1555  lastHeuristic_ = heuristic_[found]; 

1556  } 

1557  delete [] newSolution ; 

1558  newNode>setGuessedObjectiveValue(estValue) ; 

1559  tree_>push(newNode) ; 

1560  if (statistics_) { 

1561  if (numberNodes2_==maximumStatistics_) { 

1562  maximumStatistics_ = 2*maximumStatistics_; 

1563  CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_]; 

1564  memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *)); 

1565  memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *)); 

1566  delete [] statistics_; 

1567  statistics_=temp; 

1568  } 

1569  assert (!statistics_[numberNodes2_]); 

1570  statistics_[numberNodes2_]=new CbcStatistics(newNode); 

1571  } 

1572  numberNodes2_++; 

1573  # ifdef CHECK_NODE 

1574  printf("Node %x pushed on tree c\n",newNode) ; 

1575  # endif 

1576  } 

1577  else 

1578  { 

1579  if(solverCharacteristics_ && //we may be in a non standard bab 

1580  solverCharacteristics_>solutionAddsCuts()// we are in some kind of OA based bab. 

1581  ) 

1582  { 

1583  std::cerr<<"You should never get here"<<std::endl; 

1584  throw CoinError("Nodes should not be fathomed on integer infeasibility in this setting", 

1585  "branchAndBound","CbcModel") ; 

1586  } 

1587  for (i = 0 ; i < currentNumberCuts_ ; i++) 

1588  { if (addedCuts_[i]) 

1589  { if (!addedCuts_[i]>decrement(1)) 

1590  delete addedCuts_[i] ; } } 

1591  double objectiveValue = newNode>objectiveValue(); 

1592  setBestSolution(CBC_SOLUTION,objectiveValue, 

1593  solver_>getColSolution()) ; 

1594  lastHeuristic_ = NULL; 

1595  incrementUsed(solver_>getColSolution()); 

1596  //assert(nodeInfo>numberPointingToThis() <= 2) ; 

1597  // avoid accidental pruning, if newNode was final branch arm 

1598  nodeInfo>increment(); 

1599  delete newNode ; 

1600  nodeInfo>decrement() ; } } 

1601  /* 

1602  This node has been completely expanded and can be removed from the live 

1603  set. 

1604  */ 

1605  if (deleteNode) 

1606  delete node ; } 

1607  /* 

1608  End of the nonabort actions. The next block of code is executed if we've 

1609  aborted because we hit one of the limits. Clean up by deleting the live set 

1610  and break out of the node processing loop. Note that on an abort, node may 

1611  have been pushed back onto the tree for further processing, in which case 

1612  it'll be deleted in cleanTree. We need to check. 

1613  */ 

1614  else 

1615  { 

1616  tree_>cleanTree(this,COIN_DBL_MAX,bestPossibleObjective_) ; 

1617  delete nextRowCut_; 

1618  // We need to get rid of node if is has already been popped from tree 

1619  if (!nodeOnTree&&!stoppedOnGap_&&node!=rootNode) 

1620  delete node; 

1621  if (stoppedOnGap_) 

1622  { messageHandler()>message(CBC_GAP,messages()) 

1623  << bestObjective_bestPossibleObjective_ 

1624  << dblParam_[CbcAllowableGap] 

1625  << dblParam_[CbcAllowableFractionGap]*100.0 

1626  << CoinMessageEol ; 

1627  secondaryStatus_ = 2; 

1628  status_ = 0 ; } 

1629  else 

1630  if (isNodeLimitReached()) 

1631  { handler_>message(CBC_MAXNODES,messages_) << CoinMessageEol ; 

1632  secondaryStatus_ = 3; 

1633  status_ = 1 ; } 

1634  else 

1635  if (totalTime >= dblParam_[CbcMaximumSeconds]) 

1636  { handler_>message(CBC_MAXTIME,messages_) << CoinMessageEol ; 

1637  secondaryStatus_ = 4; 

1638  status_ = 1 ; } 

1639  else 

1640  if (eventHappened_) 

1641  { handler_>message(CBC_EVENT,messages_) << CoinMessageEol ; 

1642  secondaryStatus_ = 5; 

1643  status_ = 5 ; } 

1644  else 

1645  { handler_>message(CBC_MAXSOLS,messages_) << CoinMessageEol ; 

1646  secondaryStatus_ = 6; 

1647  status_ = 1 ; } 

1648  break ; } 

1649  /* 

1650  Delete cuts to get back to the original system. 

1651  

1652  I'm thinking this is redundant  the call to addCuts that conditions entry 

1653  to this code block also performs this action. 

1654  */ 

1655  int numberToDelete = getNumRows()numberRowsAtContinuous_ ; 

1656  if (numberToDelete) 

1657  { int * delRows = new int[numberToDelete] ; 

1658  int i ; 

1659  for (i = 0 ; i < numberToDelete ; i++) 

1660  { delRows[i] = i+numberRowsAtContinuous_ ; } 

1661  solver_>deleteRows(numberToDelete,delRows) ; 

1662  delete [] delRows ; } } 

1663  /* 

1664  This node fathomed when addCuts atttempted to revive it. Toss it. 

1665  */ 

1666  else 

1667  { delete node ; } } 

1668  /* 

1669  That's it, we've exhausted the search tree, or broken out of the loop because 

1670  we hit some limit on evaluation. 

1671  

1672  We may have got an intelligent tree so give it one more chance 

1673  */ 

1674  // Tell solver we are not in Branch and Cut 

1675  solver_>setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ; 

1676  tree_>endSearch(); 

1677  // If we did any sub trees  did we give up on any? 

1678  if ( numberStoppedSubTrees_) 

1679  status_=1; 

1680  if (!status_) { 

1681  bestPossibleObjective_=bestObjective_; 

1682  handler_>message(CBC_END_GOOD,messages_) 

1683  << bestObjective_ << numberIterations_ << numberNodes_<<getCurrentSeconds() 

1684  << CoinMessageEol ; 

1685  } else { 

1686  handler_>message(CBC_END,messages_) 

1687  << bestObjective_ <<bestPossibleObjective_ 

1688  << numberIterations_ << numberNodes_<<getCurrentSeconds() 

1689  << CoinMessageEol ; 

1690  } 

1691  if (numberStrongIterations_) 

1692  handler_>message(CBC_STRONG_STATS,messages_) 

1693  << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2] 

1694  << strongInfo_[1] << CoinMessageEol ; 

1695  if (statistics_) { 

1696  // report in some way 

1697  int * lookup = new int[numberObjects_]; 

1698  int i; 

1699  for (i=0;i<numberObjects_;i++) 

1700  lookup[i]=1; 

1701  bool goodIds=true; 

1702  for (i=0;i<numberObjects_;i++) { 

1703  int id = object_[i]>id(); 

1704  int iColumn = object_[i]>columnNumber(); 

1705  if (iColumn<0) 

1706  iColumn = id+numberColumns; 

1707  if(id>=0&&id<numberObjects_) { 

1708  if (lookup[id]==1) { 

1709  lookup[id]=iColumn; 

1710  } else { 

1711  goodIds=false; 

1712  break; 

1713  } 

1714  } else { 

1715  goodIds=false; 

1716  break; 

1717  } 

1718  } 

1719  if (!goodIds) { 

1720  delete [] lookup; 

1721  lookup=NULL; 

1722  } 

1723  if (doStatistics==3) { 

1724  printf(" node parent depth column value obj inf\n"); 

1725  for ( i=0;i<numberNodes2_;i++) { 

1726  statistics_[i]>print(lookup); 

1727  } 

1728  } 

1729  if (doStatistics>1) { 

1730  // Find last solution 

1731  int k; 

1732  for (k=numberNodes2_1;k>=0;k) { 

1733  if (statistics_[k]>endingObjective()!=COIN_DBL_MAX&& 

1734  !statistics_[k]>endingInfeasibility()) 

1735  break; 

1736  } 

1737  if (k>=0) { 

1738  int depth=statistics_[k]>depth(); 

1739  int * which = new int[depth+1]; 

1740  for (i=depth;i>=0;i) { 

1741  which[i]=k; 

1742  k=statistics_[k]>parentNode(); 

1743  } 

1744  printf(" node parent depth column value obj inf\n"); 

1745  for (i=0;i<=depth;i++) { 

1746  statistics_[which[i]]>print(lookup); 

1747  } 

1748  delete [] which; 

1749  } 

1750  } 

1751  // now summary 

1752  int maxDepth=0; 

1753  double averageSolutionDepth=0.0; 

1754  int numberSolutions=0; 

1755  double averageCutoffDepth=0.0; 

1756  double averageSolvedDepth=0.0; 

1757  int numberCutoff=0; 

1758  int numberDown=0; 

1759  int numberFirstDown=0; 

1760  double averageInfDown=0.0; 

1761  double averageObjDown=0.0; 

1762  int numberCutoffDown=0; 

1763  int numberUp=0; 

1764  int numberFirstUp=0; 

1765  double averageInfUp=0.0; 

1766  double averageObjUp=0.0; 

1767  int numberCutoffUp=0; 

1768  double averageNumberIterations1=0.0; 

1769  double averageValue=0.0; 

1770  for ( i=0;i<numberNodes2_;i++) { 

1771  int depth = statistics_[i]>depth(); 

1772  int way = statistics_[i]>way(); 

1773  double value = statistics_[i]>value(); 

1774  double startingObjective = statistics_[i]>startingObjective(); 

1775  int startingInfeasibility = statistics_[i]>startingInfeasibility(); 

1776  double endingObjective = statistics_[i]>endingObjective(); 

1777  int endingInfeasibility = statistics_[i]>endingInfeasibility(); 

1778  maxDepth = CoinMax(depth,maxDepth); 

1779  // Only for completed 

1780  averageNumberIterations1 += statistics_[i]>numberIterations(); 

1781  averageValue += value; 

1782  if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) { 

1783  numberSolutions++; 

1784  averageSolutionDepth += depth; 

1785  } 

1786  if (endingObjective==COIN_DBL_MAX) { 

1787  numberCutoff++; 

1788  averageCutoffDepth += depth; 

1789  if (way<0) { 

1790  numberDown++; 

1791  numberCutoffDown++; 

1792  if (way==1) 

1793  numberFirstDown++; 

1794  } else { 

1795  numberUp++; 

1796  numberCutoffUp++; 

1797  if (way==1) 

1798  numberFirstUp++; 

1799  } 

1800  } else { 

1801  averageSolvedDepth += depth; 

1802  if (way<0) { 

1803  numberDown++; 

1804  averageInfDown += startingInfeasibilityendingInfeasibility; 

1805  averageObjDown += endingObjectivestartingObjective; 

1806  if (way==1) 

1807  numberFirstDown++; 

1808  } else { 

1809  numberUp++; 

1810  averageInfUp += startingInfeasibilityendingInfeasibility; 

1811  averageObjUp += endingObjectivestartingObjective; 

1812  if (way==1) 

1813  numberFirstUp++; 

1814  } 

1815  } 

1816  } 

1817  // Now print 

1818  if (numberSolutions) 

1819  averageSolutionDepth /= (double) numberSolutions; 

1820  int numberSolved = numberNodes2_numberCutoff; 

1821  double averageNumberIterations2=numberIterations_averageNumberIterations1 

1822  numberIterationsAtContinuous; 

1823  if(numberCutoff) { 

1824  averageCutoffDepth /= (double) numberCutoff; 

1825  averageNumberIterations2 /= (double) numberCutoff; 

1826  } 

1827  if (numberNodes2_) 

1828  averageValue /= (double) numberNodes2_; 

1829  if (numberSolved) { 

1830  averageNumberIterations1 /= (double) numberSolved; 

1831  averageSolvedDepth /= (double) numberSolved; 

1832  } 

1833  printf("%d solution(s) were found (by branching) at an average depth of %g\n", 

1834  numberSolutions,averageSolutionDepth); 

1835  printf("average value of variable being branched on was %g\n", 

1836  averageValue); 

1837  printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n", 

1838  numberCutoff,averageCutoffDepth,averageNumberIterations2); 

1839  printf("%d nodes were solved at an average depth of %g with iteration count of %g\n", 

1840  numberSolved,averageSolvedDepth,averageNumberIterations1); 

1841  if (numberDown) { 

1842  averageInfDown /= (double) numberDown; 

1843  averageObjDown /= (double) numberDown; 

1844  } 

1845  printf("Down %d nodes (%d first, %d second)  %d cutoff, rest decrease numinf %g increase obj %g\n", 

1846  numberDown,numberFirstDown,numberDownnumberFirstDown,numberCutoffDown, 

1847  averageInfDown,averageObjDown); 

1848  if (numberUp) { 

1849  averageInfUp /= (double) numberUp; 

1850  averageObjUp /= (double) numberUp; 

1851  } 

1852  printf("Up %d nodes (%d first, %d second)  %d cutoff, rest decrease numinf %g increase obj %g\n", 

1853  numberUp,numberFirstUp,numberUpnumberFirstUp,numberCutoffUp, 

1854  averageInfUp,averageObjUp); 

1855  for ( i=0;i<numberNodes2_;i++) 

1856  delete statistics_[i]; 

1857  delete [] statistics_; 

1858  statistics_=NULL; 

1859  maximumStatistics_=0; 

1860  delete [] lookup; 

1861  } 

1862  /* 

1863  If we think we have a solution, restore and confirm it with a call to 

1864  setBestSolution(). We need to reset the cutoff value so as not to fathom 

1865  the solution on bounds. Note that calling setBestSolution( ..., true) 

1866  leaves the continuousSolver_ bounds vectors fixed at the solution value. 

1867  

1868  Running resolve() here is a failsafe  setBestSolution has already 

1869  reoptimised using the continuousSolver_. If for some reason we fail to 

1870  prove optimality, run the problem again after instructing the solver to 

1871  tell us more. 

1872  

1873  If all looks good, replace solver_ with continuousSolver_, so that the 

1874  outside world will be able to obtain information about the solution using 

1875  public methods. 

1876  */ 

1877  if (bestSolution_&&solverCharacteristics_>solverType()<2) 

1878  { setCutoff(1.0e50) ; // As best solution should be worse than cutoff 

1879  phase_=5; 

1880  double increment = getDblParam(CbcModel::CbcCutoffIncrement) ; 

1881  bestObjective_ += 100.0*increment+1.0e3; 

1882  setBestSolution(CBC_END_SOLUTION,bestObjective_,bestSolution_,true) ; 

1883  continuousSolver_>resolve() ; 

1884  if (!continuousSolver_>isProvenOptimal()) 

1885  { continuousSolver_>messageHandler()>setLogLevel(2) ; 

1886  continuousSolver_>initialSolve() ; } 

1887  delete solver_ ; 

1888  // above deletes solverCharacteristics_ 

1889  solverCharacteristics_ = NULL; 

1890  solver_ = continuousSolver_ ; 

1891  setPointers(solver_); 

1892  continuousSolver_ = NULL ; } 

1893  /* 

1894  Clean up dangling objects. continuousSolver_ may already be toast. 

1895  */ 

1896  delete lastws ; 

1897  delete [] whichGenerator_ ; 

1898  whichGenerator_=NULL; 

1899  delete [] lowerBefore ; 

1900  delete [] upperBefore ; 

1901  delete [] walkback_ ; 

1902  walkback_ = NULL ; 

1903  delete [] addedCuts_ ; 

1904  addedCuts_ = NULL ; 

1905  // Get rid of characteristics 

1906  solverCharacteristics_=NULL; 

1907  if (continuousSolver_) 

1908  { delete continuousSolver_ ; 

1909  continuousSolver_ = NULL ; } 

1910  /* 

1911  Destroy global cuts by replacing with an empty OsiCuts object. 

1912  */ 

1913  globalCuts_= OsiCuts() ; 

1914  if (strategy_&&strategy_>preProcessState()>0) { 

1915  // undo preprocessing 

1916  CglPreProcess * process = strategy_>process(); 

1917  assert (process); 

1918  int n = originalSolver>getNumCols(); 

1919  if (bestSolution_) { 

1920  delete [] bestSolution_; 

1921  bestSolution_ = new double [n]; 

1922  process>postProcess(*solver_); 

1923  } 

1924  strategy_>deletePreProcess(); 

1925  // Solution now back in originalSolver 

1926  delete solver_; 

1927  solver_=originalSolver; 

1928  if (bestSolution_) 

1929  memcpy(bestSolution_,solver_>getColSolution(),n*sizeof(double)); 

1930  // put back original objects if there were any 

1931  if (originalObject) { 

1932  int iColumn; 

1933  for (iColumn=0;iColumn<numberObjects_;iColumn++) 

1934  delete object_[iColumn]; 

1935  delete [] object_; 

1936  numberObjects_ = numberOriginalObjects; 

1937  object_=originalObject; 

1938  delete [] integerVariable_; 

1939  numberIntegers_=0; 

1940  for (iColumn=0;iColumn<n;iColumn++) { 

1941  if (solver_>isInteger(iColumn)) 

1942  numberIntegers_++; 

1943  } 

1944  integerVariable_ = new int[numberIntegers_]; 

1945  numberIntegers_=0; 

1946  for (iColumn=0;iColumn<n;iColumn++) { 

1947  if (solver_>isInteger(iColumn)) 

1948  integerVariable_[numberIntegers_++]=iColumn; 

1949  } 

1950  } 

1951  } 

1952  return ; } 

1953  

1954  

1955  // Solve the initial LP relaxation 

1956  void 

1957  CbcModel::initialSolve() 

1958  { 

1959  assert (solver_); 

1960  assert (!solverCharacteristics_); 

1961  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_>getAuxiliaryInfo()); 

1962  if (solverCharacteristics) { 

1963  solverCharacteristics_ = solverCharacteristics; 

1964  } else { 

1965  // replace in solver 

1966  OsiBabSolver defaultC; 

1967  solver_>setAuxiliaryInfo(&defaultC); 

1968  solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_>getAuxiliaryInfo()); 

1969  } 

1970  solverCharacteristics_>setSolver(solver_); 

1971  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ; 

1972  solver_>initialSolve(); 

1973  solver_>setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ; 

1974  // But set up so Jon Lee will be happy 

1975  status_=1; 

1976  secondaryStatus_ = 1; 

1977  originalContinuousObjective_ = solver_>getObjValue()*solver_>getObjSense(); 

1978  delete [] continuousSolution_; 

1979  continuousSolution_ = CoinCopyOfArray(solver_>getColSolution(), 

1980  solver_>getNumCols()); 

1981  setPointers(solver_); 

1982  solverCharacteristics_ = NULL; 

1983  } 

1984  

1985  /*! \brief Get an empty basis object 

1986  

1987  Return an empty CoinWarmStartBasis object with the requested capacity, 

1988  appropriate for the current solver. The object is cloned from the object 

1989  cached as emptyWarmStart_. If there is no cached object, the routine 

1990  queries the solver for a warm start object, empties it, and caches the 

1991  result. 

1992  */ 

1993  

1994  CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const 

1995  

1996  { CoinWarmStartBasis *emptyBasis ; 

1997  /* 

1998  Acquire an empty basis object, if we don't yet have one. 

1999  */ 

2000  if (emptyWarmStart_ == 0) 

2001  { if (solver_ == 0) 

2002  { throw CoinError("Cannot construct basis without solver!", 

2003  "getEmptyBasis","CbcModel") ; } 

2004  emptyBasis = 

2005  dynamic_cast<CoinWarmStartBasis *>(solver_>getEmptyWarmStart()) ; 

2006  if (emptyBasis == 0) 

2007  { throw CoinError( 

2008  "Solver does not appear to use a basisoriented warm start.", 

2009  "getEmptyBasis","CbcModel") ; } 

2010  emptyBasis>setSize(0,0) ; 

2011  emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ; } 

2012  /* 

2013  Clone the empty basis object, resize it as requested, and return. 

2014  */ 

2015  emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_>clone()) ; 

2016  assert(emptyBasis) ; 

2017  if (ns != 0  na != 0) emptyBasis>setSize(ns,na) ; 

2018  

2019  return (emptyBasis) ; } 

2020  

2021  

2022  /** Default Constructor 

2023  

2024  Creates an empty model without an associated solver. 

2025  */ 

2026  CbcModel::CbcModel() 

2027  

2028  : 

2029  solver_(NULL), 

2030  ourSolver_(true), 

2031  continuousSolver_(NULL), 

2032  referenceSolver_(NULL), 

2033  defaultHandler_(true), 

2034  emptyWarmStart_(NULL), 

2035  bestObjective_(COIN_DBL_MAX), 

2036  bestPossibleObjective_(COIN_DBL_MAX), 

2037  sumChangeObjective1_(0.0), 

2038  sumChangeObjective2_(0.0), 

2039  bestSolution_(NULL), 

2040  currentSolution_(NULL), 

2041  testSolution_(NULL), 

2042  minimumDrop_(1.0e4), 

2043  numberSolutions_(0), 

2044  stateOfSearch_(0), 

2045  hotstartSolution_(NULL), 

2046  hotstartPriorities_(NULL), 

2047  numberHeuristicSolutions_(0), 

2048  numberNodes_(0), 

2049  numberNodes2_(0), 

2050  numberIterations_(0), 

2051  status_(1), 

2052  secondaryStatus_(1), 

2053  numberIntegers_(0), 

2054  numberRowsAtContinuous_(0), 

2055  maximumNumberCuts_(0), 

2056  phase_(0), 

2057  currentNumberCuts_(0), 

2058  maximumDepth_(0), 

2059  walkback_(NULL), 

2060  addedCuts_(NULL), 

2061  nextRowCut_(NULL), 

2062  currentNode_(NULL), 

2063  integerVariable_(NULL), 

2064  integerInfo_(NULL), 

2065  continuousSolution_(NULL), 

2066  usedInSolution_(NULL), 

2067  specialOptions_(0), 

2068  subTreeModel_(NULL), 

2069  numberStoppedSubTrees_(0), 

2070  presolve_(0), 

2071  numberStrong_(5), 

2072  numberBeforeTrust_(0), 

2073  numberPenalties_(20), 

2074  penaltyScaleFactor_(3.0), 

2075  numberAnalyzeIterations_(0), 

2076  analyzeResults_(NULL), 

2077  numberInfeasibleNodes_(0), 

2078  problemType_(0), 

2079  printFrequency_(0), 

2080  numberCutGenerators_(0), 

2081  generator_(NULL), 

2082  virginGenerator_(NULL), 

2083  numberHeuristics_(0), 

2084  heuristic_(NULL), 

2085  lastHeuristic_(NULL), 

2086  eventHandler_(0), 

2087  numberObjects_(0), 

2088  object_(NULL), 

2089  originalColumns_(NULL), 

2090  howOftenGlobalScan_(1), 

2091  numberGlobalViolations_(0), 

2092  continuousObjective_(COIN_DBL_MAX), 

2093  originalContinuousObjective_(COIN_DBL_MAX), 

2094  continuousInfeasibilities_(INT_MAX), 

2095  maximumCutPassesAtRoot_(20), 

2096  maximumCutPasses_(10), 

2097  currentPassNumber_(0), 

2098  maximumWhich_(1000), 

2099  whichGenerator_(NULL), 

2100  maximumStatistics_(0), 

2101  statistics_(NULL), 

2102  numberFixedAtRoot_(0), 

2103  numberFixedNow_(0), 

2104  stoppedOnGap_(false), 

2105  eventHappened_(false), 

2106  numberLongStrong_(0), 

2107  numberOldActiveCuts_(0), 

2108  numberNewCuts_(0), 

2109  sizeMiniTree_(0), 

2110  searchStrategy_(1), 

2111  numberStrongIterations_(0), 

2112  resolveAfterTakeOffCuts_(true) 

2113  { 

2114  intParam_[CbcMaxNumNode] = 2147483647; 

2115  intParam_[CbcMaxNumSol] = 9999999; 

2116  intParam_[CbcFathomDiscipline] = 0; 

2117  

2118  dblParam_[CbcIntegerTolerance] = 1e6; 

2119  dblParam_[CbcInfeasibilityWeight] = 0.0; 

2120  dblParam_[CbcCutoffIncrement] = 1e5; 

2121  dblParam_[CbcAllowableGap] = 1.0e10; 

2122  dblParam_[CbcAllowableFractionGap] = 0.0; 

2123  dblParam_[CbcMaximumSeconds] = 1.0e100; 

2124  dblParam_[CbcCurrentCutoff] = 1.0e100; 

2125  dblParam_[CbcOptimizationDirection] = 1.0; 

2126  dblParam_[CbcCurrentObjectiveValue] = 1.0e100; 

2127  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100; 

2128  dblParam_[CbcStartSeconds] = 0.0; 

2129  strongInfo_[0]=0; 

2130  strongInfo_[1]=0; 

2131  strongInfo_[2]=0; 

2132  solverCharacteristics_ = NULL; 

2133  nodeCompare_=new CbcCompareDefault();; 

2134  problemFeasibility_=new CbcFeasibilityBase(); 

2135  tree_= new CbcTree(); 

2136  branchingMethod_=NULL; 

2137  strategy_=NULL; 

2138  parentModel_=NULL; 

2139  cbcColLower_ = NULL; 

2140  cbcColUpper_ = NULL; 

2141  cbcRowLower_ = NULL; 

2142  cbcRowUpper_ = NULL; 

2143  cbcColSolution_ = NULL; 

2144  cbcRowPrice_ = NULL; 

2145  cbcReducedCost_ = NULL; 

2146  cbcRowActivity_ = NULL; 

2147  appData_=NULL; 

2148  handler_ = new CoinMessageHandler(); 

2149  handler_>setLogLevel(2); 

2150  messages_ = CbcMessage(); 

2151  eventHandler_ = new CbcEventHandler() ; 

2152  } 

2153  

2154  /** Constructor from solver. 

2155  

2156  Creates a model complete with a clone of the solver passed as a parameter. 

2157  */ 

2158  

2159  CbcModel::CbcModel(const OsiSolverInterface &rhs) 

2160  : 

2161  continuousSolver_(NULL), 

2162  referenceSolver_(NULL), 

2163  defaultHandler_(true), 

2164  emptyWarmStart_(NULL), 

2165  bestObjective_(COIN_DBL_MAX), 

2166  bestPossibleObjective_(COIN_DBL_MAX), 

2167  sumChangeObjective1_(0.0), 

2168  sumChangeObjective2_(0.0), 

2169  minimumDrop_(1.0e4), 

2170  numberSolutions_(0), 

2171  stateOfSearch_(0), 

2172  hotstartSolution_(NULL), 

2173  hotstartPriorities_(NULL), 

2174  numberHeuristicSolutions_(0), 

2175  numberNodes_(0), 

2176  numberNodes2_(0), 

2177  numberIterations_(0), 

2178  status_(1), 

2179  secondaryStatus_(1), 

2180  numberRowsAtContinuous_(0), 

2181  maximumNumberCuts_(0), 

2182  phase_(0), 

2183  currentNumberCuts_(0), 

2184  maximumDepth_(0), 

2185  walkback_(NULL), 

2186  addedCuts_(NULL), 

2187  nextRowCut_(NULL), 

2188  currentNode_(NULL), 

2189  integerInfo_(NULL), 

2190  specialOptions_(0), 

2191  subTreeModel_(NULL), 

2192  numberStoppedSubTrees_(0), 

2193  presolve_(0), 

2194  numberStrong_(5), 

2195  numberBeforeTrust_(0), 

2196  numberPenalties_(20), 

2197  penaltyScaleFactor_(3.0), 

2198  numberAnalyzeIterations_(0), 

2199  analyzeResults_(NULL), 

2200  numberInfeasibleNodes_(0), 

2201  problemType_(0), 

2202  printFrequency_(0), 

2203  numberCutGenerators_(0), 

2204  generator_(NULL), 

2205  virginGenerator_(NULL), 

2206  numberHeuristics_(0), 

2207  heuristic_(NULL), 

2208  lastHeuristic_(NULL), 

2209  eventHandler_(0), 

2210  numberObjects_(0), 

2211  object_(NULL), 

2212  originalColumns_(NULL), 

2213  howOftenGlobalScan_(1), 

2214  numberGlobalViolations_(0), 

2215  continuousObjective_(COIN_DBL_MAX), 

2216  originalContinuousObjective_(COIN_DBL_MAX), 

2217  continuousInfeasibilities_(INT_MAX), 

2218  maximumCutPassesAtRoot_(20), 

2219  maximumCutPasses_(10), 

2220  currentPassNumber_(0), 

2221  maximumWhich_(1000), 

2222  whichGenerator_(NULL), 

2223  maximumStatistics_(0), 

2224  statistics_(NULL), 

2225  numberFixedAtRoot_(0), 

2226  numberFixedNow_(0), 

2227  stoppedOnGap_(false), 

2228  eventHappened_(false), 

2229  numberLongStrong_(0), 

2230  numberOldActiveCuts_(0), 

2231  numberNewCuts_(0), 

2232  sizeMiniTree_(0), 

2233  searchStrategy_(1), 

2234  numberStrongIterations_(0), 

2235  resolveAfterTakeOffCuts_(true) 

2236  { 

2237  intParam_[CbcMaxNumNode] = 2147483647; 

2238  intParam_[CbcMaxNumSol] = 9999999; 

2239  intParam_[CbcFathomDiscipline] = 0; 

2240  

2241  dblParam_[CbcIntegerTolerance] = 1e6; 

2242  dblParam_[CbcInfeasibilityWeight] = 0.0; 

2243  dblParam_[CbcCutoffIncrement] = 1e5; 

2244  dblParam_[CbcAllowableGap] = 1.0e10; 

2245  dblParam_[CbcAllowableFractionGap] = 0.0; 

2246  dblParam_[CbcMaximumSeconds] = 1.0e100; 

2247  dblParam_[CbcCurrentCutoff] = 1.0e100; 

2248  dblParam_[CbcOptimizationDirection] = 1.0; 

2249  dblParam_[CbcCurrentObjectiveValue] = 1.0e100; 

2250  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100; 

2251  dblParam_[CbcStartSeconds] = 0.0; 

2252  strongInfo_[0]=0; 

2253  strongInfo_[1]=0; 

2254  strongInfo_[2]=0; 

2255  solverCharacteristics_ = NULL; 

2256  

2257  nodeCompare_=new CbcCompareDefault();; 

2258  problemFeasibility_=new CbcFeasibilityBase(); 

2259  tree_= new CbcTree(); 

2260  branchingMethod_=NULL; 

2261  strategy_=NULL; 

2262  parentModel_=NULL; 

2263  appData_=NULL; 

2264  handler_ = new CoinMessageHandler(); 

2265  handler_>setLogLevel(2); 

2266  messages_ = CbcMessage(); 

2267  eventHandler_ = new CbcEventHandler() ; 

2268  solver_ = rhs.clone(); 

2269  referenceSolver_ = solver_>clone(); 

2270  ourSolver_ = true ; 

2271  cbcColLower_ = NULL; 

2272  cbcColUpper_ = NULL; 

2273  cbcRowLower_ = NULL; 

2274  cbcRowUpper_ = NULL; 

2275  cbcColSolution_ = NULL; 

2276  cbcRowPrice_ = NULL; 

2277  cbcReducedCost_ = NULL; 

2278  cbcRowActivity_ = NULL; 

2279  

2280  // Initialize solution and integer variable vectors 

2281  bestSolution_ = NULL; // to say no solution found 

2282  numberIntegers_=0; 

2283  int numberColumns = solver_>getNumCols(); 

2284  int iColumn; 

2285  if (numberColumns) { 

2286  // Space for current solution 

2287  currentSolution_ = new double[numberColumns]; 

2288  continuousSolution_ = new double[numberColumns]; 

2289  usedInSolution_ = new int[numberColumns]; 

2290  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

2291  if( solver_>isInteger(iColumn)) 

2292  numberIntegers_++; 

2293  } 

2294  } else { 

2295  // empty model 

2296  currentSolution_=NULL; 

2297  continuousSolution_=NULL; 

2298  usedInSolution_=NULL; 

2299  } 

2300  testSolution_=currentSolution_; 

2301  if (numberIntegers_) { 

2302  integerVariable_ = new int [numberIntegers_]; 

2303  numberIntegers_=0; 

2304  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

2305  if( solver_>isInteger(iColumn)) 

2306  integerVariable_[numberIntegers_++]=iColumn; 

2307  } 

2308  } else { 

2309  integerVariable_ = NULL; 

2310  } 

2311  } 

2312  

2313  /* 

2314  Assign a solver to the model (model assumes ownership) 

2315  

2316  The integer variable vector is initialized if it's not already present. 

2317  If deleteSolver then current solver deleted (if model owned) 

2318  

2319  Assuming ownership matches usage in OsiSolverInterface 

2320  (cf. assignProblem, loadProblem). 

2321  

2322  TODO: What to do about solver parameters? A simple copy likely won't do it, 

2323  because the SI must push the settings into the underlying solver. In 

2324  the context of switching solvers in cbc, this means that command line 

2325  settings will get lost. Stash the command line somewhere and reread it 

2326  here, maybe? 

2327  

2328  TODO: More generally, how much state should be transferred from the old 

2329  solver to the new solver? Best perhaps to see how usage develops. 

2330  What's done here mimics the CbcModel(OsiSolverInterface) constructor. 

2331  */ 

2332  void 

2333  CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver) 

2334  

2335  { 

2336  // Keep the current message level for solver (if solver exists) 

2337  if (solver_) 

2338  solver>messageHandler()>setLogLevel(solver_>messageHandler()>logLevel()) ; 

2339  

2340  if (ourSolver_&&deleteSolver) delete solver_ ; 

2341  solver_ = solver; 

2342  solver = NULL ; 

2343  ourSolver_ = true ; 

2344  /* 

2345  Basis information is solverspecific. 

2346  */ 

2347  if (emptyWarmStart_) 

2348  { delete emptyWarmStart_ ; 

2349  emptyWarmStart_ = 0 ; } 

2350  /* 

2351  Initialize integer variable vector. 

2352  */ 

2353  numberIntegers_=0; 

2354  int numberColumns = solver_>getNumCols(); 

2355  int iColumn; 

2356  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

2357  if( solver_>isInteger(iColumn)) 

2358  numberIntegers_++; 

2359  } 

2360  if (numberIntegers_) { 

2361  delete [] integerVariable_; 

2362  integerVariable_ = new int [numberIntegers_]; 

2363  numberIntegers_=0; 

2364  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

2365  if( solver_>isInteger(iColumn)) 

2366  integerVariable_[numberIntegers_++]=iColumn; 

2367  } 

2368  } else { 

2369  integerVariable_ = NULL; 

2370  } 

2371  

2372  return ; 

2373  } 

2374  

2375  // Copy constructor. 

2376  

2377  CbcModel::CbcModel(const CbcModel & rhs, bool noTree) 

2378  : 

2379  continuousSolver_(NULL), 

2380  referenceSolver_(NULL), 

2381  defaultHandler_(rhs.defaultHandler_), 

2382  emptyWarmStart_(NULL), 

2383  bestObjective_(rhs.bestObjective_), 

2384  bestPossibleObjective_(rhs.bestPossibleObjective_), 

2385  sumChangeObjective1_(rhs.sumChangeObjective1_), 

2386  sumChangeObjective2_(rhs.sumChangeObjective2_), 

2387  minimumDrop_(rhs.minimumDrop_), 

2388  numberSolutions_(rhs.numberSolutions_), 

2389  stateOfSearch_(rhs.stateOfSearch_), 

2390  numberHeuristicSolutions_(rhs.numberHeuristicSolutions_), 

2391  numberNodes_(rhs.numberNodes_), 

2392  numberNodes2_(rhs.numberNodes2_), 

2393  numberIterations_(rhs.numberIterations_), 

2394  status_(rhs.status_), 

2395  secondaryStatus_(rhs.secondaryStatus_), 

2396  specialOptions_(rhs.specialOptions_), 

2397  subTreeModel_(rhs.subTreeModel_), 

2398  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_), 

2399  presolve_(rhs.presolve_), 

2400  numberStrong_(rhs.numberStrong_), 

2401  numberBeforeTrust_(rhs.numberBeforeTrust_), 

2402  numberPenalties_(rhs.numberPenalties_), 

2403  penaltyScaleFactor_(rhs.penaltyScaleFactor_), 

2404  numberAnalyzeIterations_(rhs.numberAnalyzeIterations_), 

2405  analyzeResults_(NULL), 

2406  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_), 

2407  problemType_(rhs.problemType_), 

2408  printFrequency_(rhs.printFrequency_), 

2409  howOftenGlobalScan_(rhs.howOftenGlobalScan_), 

2410  numberGlobalViolations_(rhs.numberGlobalViolations_), 

2411  continuousObjective_(rhs.continuousObjective_), 

2412  originalContinuousObjective_(rhs.originalContinuousObjective_), 

2413  continuousInfeasibilities_(rhs.continuousInfeasibilities_), 

2414  maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_), 

2415  maximumCutPasses_( rhs.maximumCutPasses_), 

2416  currentPassNumber_(rhs.currentPassNumber_), 

2417  maximumWhich_(rhs.maximumWhich_), 

2418  whichGenerator_(NULL), 

2419  maximumStatistics_(0), 

2420  statistics_(NULL), 

2421  numberFixedAtRoot_(rhs.numberFixedAtRoot_), 

2422  numberFixedNow_(rhs.numberFixedNow_), 

2423  stoppedOnGap_(rhs.stoppedOnGap_), 

2424  eventHappened_(rhs.eventHappened_), 

2425  numberLongStrong_(rhs.numberLongStrong_), 

2426  numberOldActiveCuts_(rhs.numberOldActiveCuts_), 

2427  numberNewCuts_(rhs.numberNewCuts_), 

2428  sizeMiniTree_(rhs.sizeMiniTree_), 

2429  searchStrategy_(rhs.searchStrategy_), 

2430  numberStrongIterations_(rhs.numberStrongIterations_), 

2431  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_) 

2432  { 

2433  intParam_[CbcMaxNumNode] = rhs.intParam_[CbcMaxNumNode]; 

2434  intParam_[CbcMaxNumSol] = rhs.intParam_[CbcMaxNumSol]; 

2435  intParam_[CbcFathomDiscipline] = rhs.intParam_[CbcFathomDiscipline]; 

2436  dblParam_[CbcIntegerTolerance] = rhs.dblParam_[CbcIntegerTolerance]; 

2437  dblParam_[CbcInfeasibilityWeight] = rhs.dblParam_[CbcInfeasibilityWeight]; 

2438  dblParam_[CbcCutoffIncrement] = rhs.dblParam_[CbcCutoffIncrement]; 

2439  dblParam_[CbcAllowableGap] = rhs.dblParam_[CbcAllowableGap]; 

2440  dblParam_[CbcAllowableFractionGap] = rhs.dblParam_[CbcAllowableFractionGap]; 

2441  dblParam_[CbcMaximumSeconds] = rhs.dblParam_[CbcMaximumSeconds]; 

2442  dblParam_[CbcCurrentCutoff] = rhs.dblParam_[CbcCurrentCutoff]; 

2443  dblParam_[CbcOptimizationDirection] = rhs.dblParam_[CbcOptimizationDirection]; 

2444  dblParam_[CbcCurrentObjectiveValue] = rhs.dblParam_[CbcCurrentObjectiveValue]; 

2445  dblParam_[CbcCurrentMinimizationObjectiveValue] = rhs.dblParam_[CbcCurrentMinimizationObjectiveValue]; 

2446  dblParam_[CbcStartSeconds] = dblParam_[CbcStartSeconds]; // will be overwritten hopefully 

2447  strongInfo_[0]=rhs.strongInfo_[0]; 

2448  strongInfo_[1]=rhs.strongInfo_[1]; 

2449  strongInfo_[2]=rhs.strongInfo_[2]; 

2450  solverCharacteristics_ = NULL; 

2451  if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_>clone() ; 

2452  if (defaultHandler_) { 

2453  handler_ = new CoinMessageHandler(); 

2454  handler_>setLogLevel(2); 

2455  } else { 

2456  handler_ = rhs.handler_; 

2457  } 

2458  messageHandler()>setLogLevel(rhs.messageHandler()>logLevel()); 

2459  numberCutGenerators_ = rhs.numberCutGenerators_; 

2460  if (numberCutGenerators_) { 

2461  generator_ = new CbcCutGenerator * [numberCutGenerators_]; 

2462  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_]; 

2463  int i; 

2464  for (i=0;i<numberCutGenerators_;i++) { 

2465  generator_[i]=new CbcCutGenerator(*rhs.generator_[i]); 

2466  virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]); 

2467  } 

2468  } else { 

2469  generator_=NULL; 

2470  virginGenerator_=NULL; 

2471  } 

2472  if (!noTree) 

2473  globalCuts_ = rhs.globalCuts_; 

2474  numberHeuristics_ = rhs.numberHeuristics_; 

2475  if (numberHeuristics_) { 

2476  heuristic_ = new CbcHeuristic * [numberHeuristics_]; 

2477  int i; 

2478  for (i=0;i<numberHeuristics_;i++) { 

2479  heuristic_[i]=rhs.heuristic_[i]>clone(); 

2480  } 

2481  } else { 

2482  heuristic_=NULL; 

2483  } 

2484  lastHeuristic_ = NULL; 

2485  if (rhs.eventHandler_) 

2486  { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; } 

2487  else 

2488  { eventHandler_ = NULL ; } 

2489  numberObjects_=rhs.numberObjects_; 

2490  if (numberObjects_) { 

2491  object_ = new CbcObject * [numberObjects_]; 

2492  int i; 

2493  for (i=0;i<numberObjects_;i++) 

2494  object_[i]=(rhs.object_[i])>clone(); 

2495  } else { 

2496  object_=NULL; 

2497  } 

2498  if (rhs.referenceSolver_) 

2499  referenceSolver_ = rhs.referenceSolver_>clone(); 

2500  else 

2501  referenceSolver_=NULL; 

2502  if (!noTree!rhs.continuousSolver_) 

2503  solver_ = rhs.solver_>clone(); 

2504  else 

2505  solver_ = rhs.continuousSolver_>clone(); 

2506  if (rhs.originalColumns_) { 

2507  int numberColumns = solver_>getNumCols(); 

2508  originalColumns_= new int [numberColumns]; 

2509  memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int)); 

2510  } else { 

2511  originalColumns_=NULL; 

2512  } 

2513  if (maximumWhich_&&rhs.whichGenerator_) 

2514  whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_); 

2515  nodeCompare_=rhs.nodeCompare_>clone(); 

2516  problemFeasibility_=rhs.problemFeasibility_>clone(); 

2517  tree_= rhs.tree_>clone(); 

2518  branchingMethod_=rhs.branchingMethod_; 

2519  cbcColLower_ = NULL; 

2520  cbcColUpper_ = NULL; 

2521  cbcRowLower_ = NULL; 

2522  cbcRowUpper_ = NULL; 

2523  cbcColSolution_ = NULL; 

2524  cbcRowPrice_ = NULL; 

2525  cbcReducedCost_ = NULL; 

2526  cbcRowActivity_ = NULL; 

2527  if (rhs.strategy_) 

2528  strategy_=rhs.strategy_>clone(); 

2529  else 

2530  strategy_=NULL; 

2531  parentModel_=rhs.parentModel_; 

2532  appData_=rhs.appData_; 

2533  messages_ = rhs.messages_; 

2534  ourSolver_ = true ; 

2535  messageHandler()>setLogLevel(rhs.messageHandler()>logLevel()); 

2536  numberIntegers_=rhs.numberIntegers_; 

2537  if (numberIntegers_) { 

2538  integerVariable_ = new int [numberIntegers_]; 

2539  memcpy(integerVariable_,rhs.integerVariable_,numberIntegers_*sizeof(int)); 

2540  integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_>getNumCols()); 

2541  } else { 

2542  integerVariable_ = NULL; 

2543  integerInfo_=NULL; 

2544  } 

2545  if (rhs.hotstartSolution_) { 

2546  int numberColumns = solver_>getNumCols(); 

2547  hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns); 

2548  hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns); 

2549  } else { 

2550  hotstartSolution_ = NULL; 

2551  hotstartPriorities_ =NULL; 

2552  } 

2553  if (rhs.bestSolution_&&!noTree) { 

2554  int numberColumns = solver_>getNumCols(); 

2555  bestSolution_ = new double[numberColumns]; 

2556  memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double)); 

2557  } else { 

2558  bestSolution_=NULL; 

2559  } 

2560  if (!noTree) { 

2561  int numberColumns = solver_>getNumCols(); 

2562  currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns); 

2563  continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns); 

2564  usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns); 

2565  } else { 

2566  currentSolution_=NULL; 

2567  continuousSolution_=NULL; 

2568  usedInSolution_=NULL; 

2569  } 

2570  testSolution_=currentSolution_; 

2571  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_; 

2572  maximumNumberCuts_=rhs.maximumNumberCuts_; 

2573  phase_ = rhs.phase_; 

2574  currentNumberCuts_=rhs.currentNumberCuts_; 

2575  maximumDepth_= rhs.maximumDepth_; 

2576  if (noTree) { 

2577  bestObjective_ = COIN_DBL_MAX; 

2578  numberSolutions_ =0; 

2579  stateOfSearch_= 0; 

2580  numberHeuristicSolutions_=0; 

2581  numberNodes_=0; 

2582  numberNodes2_=0; 

2583  numberIterations_=0; 

2584  status_=0; 

2585  subTreeModel_=NULL; 

2586  numberStoppedSubTrees_=0; 

2587  continuousObjective_=COIN_DBL_MAX; 

2588  originalContinuousObjective_=COIN_DBL_MAX; 

2589  continuousInfeasibilities_=INT_MAX; 

2590  maximumNumberCuts_=0; 

2591  tree_>cleanTree(this,COIN_DBL_MAX,bestPossibleObjective_); 

2592  bestPossibleObjective_ = COIN_DBL_MAX; 

2593  } 

2594  // These are only used as temporary arrays so need not be filled 

2595  if (maximumNumberCuts_) { 

2596  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

2597  } else { 

2598  addedCuts_ = NULL; 

2599  } 

2600  nextRowCut_ = NULL; 

2601  currentNode_ = NULL; 

2602  if (maximumDepth_) 

2603  walkback_ = new CbcNodeInfo * [maximumDepth_]; 

2604  else 

2605  walkback_ = NULL; 

2606  synchronizeModel(); 

2607  } 

2608  

2609  // Assignment operator 

2610  CbcModel & 

2611  CbcModel::operator=(const CbcModel& rhs) 

2612  { 

2613  if (this!=&rhs) { 

2614  gutsOfDestructor(); 

2615  if (defaultHandler_) 

2616  { delete handler_; 

2617  handler_ = NULL; } 

2618  defaultHandler_ = rhs.defaultHandler_; 

2619  if (defaultHandler_) 

2620  { handler_ = new CoinMessageHandler(); 

2621  handler_>setLogLevel(2); } 

2622  else 

2623  { handler_ = rhs.handler_; } 

2624  messages_ = rhs.messages_; 

2625  messageHandler()>setLogLevel(rhs.messageHandler()>logLevel()); 

2626  

2627  delete solver_; 

2628  if (rhs.solver_) 

2629  { solver_ = rhs.solver_>clone() ; } 

2630  else 

2631  { solver_ = 0 ; } 

2632  ourSolver_ = true ; 

2633  delete continuousSolver_ ; 

2634  if (rhs.continuousSolver_) 

2635  { solver_ = rhs.continuousSolver_>clone() ; } 

2636  else 

2637  { continuousSolver_ = 0 ; } 

2638  

2639  if (rhs.referenceSolver_) 

2640  { solver_ = rhs.referenceSolver_>clone() ; } 

2641  else 

2642  { referenceSolver_ = NULL ; } 

2643  

2644  delete emptyWarmStart_ ; 

2645  if (rhs.emptyWarmStart_) 

2646  { emptyWarmStart_ = rhs.emptyWarmStart_>clone() ; } 

2647  else 

2648  { emptyWarmStart_ = 0 ; } 

2649  

2650  bestObjective_ = rhs.bestObjective_; 

2651  bestPossibleObjective_=rhs.bestPossibleObjective_; 

2652  sumChangeObjective1_=rhs.sumChangeObjective1_; 

2653  sumChangeObjective2_=rhs.sumChangeObjective2_; 

2654  delete [] bestSolution_; 

2655  if (rhs.bestSolution_) { 

2656  int numberColumns = rhs.getNumCols(); 

2657  bestSolution_ = new double[numberColumns]; 

2658  memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double)); 

2659  } else { 

2660  bestSolution_=NULL; 

2661  } 

2662  int numberColumns = solver_>getNumCols(); 

2663  currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns); 

2664  continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns); 

2665  usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns); 

2666  testSolution_=currentSolution_; 

2667  minimumDrop_ = rhs.minimumDrop_; 

2668  numberSolutions_=rhs.numberSolutions_; 

2669  stateOfSearch_= rhs.stateOfSearch_; 

2670  numberHeuristicSolutions_=rhs.numberHeuristicSolutions_; 

2671  numberNodes_ = rhs.numberNodes_; 

2672  numberNodes2_ = rhs.numberNodes2_; 

2673  numberIterations_ = rhs.numberIterations_; 

2674  status_ = rhs.status_; 

2675  secondaryStatus_ = rhs.secondaryStatus_; 

2676  specialOptions_ = rhs.specialOptions_; 

2677  subTreeModel_ = rhs.subTreeModel_; 

2678  numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_; 

2679  presolve_ = rhs.presolve_; 

2680  numberStrong_ = rhs.numberStrong_; 

2681  numberBeforeTrust_ = rhs.numberBeforeTrust_; 

2682  numberPenalties_ = rhs.numberPenalties_; 

2683  penaltyScaleFactor_ = rhs.penaltyScaleFactor_; 

2684  numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_; 

2685  delete [] analyzeResults_; 

2686  analyzeResults_ = NULL; 

2687  numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_; 

2688  problemType_ = rhs.problemType_; 

2689  printFrequency_ = rhs.printFrequency_; 

2690  howOftenGlobalScan_=rhs.howOftenGlobalScan_; 

2691  numberGlobalViolations_=rhs.numberGlobalViolations_; 

2692  continuousObjective_=rhs.continuousObjective_; 

2693  originalContinuousObjective_ = rhs.originalContinuousObjective_; 

2694  continuousInfeasibilities_ = rhs.continuousInfeasibilities_; 

2695  maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_; 

2696  maximumCutPasses_ = rhs.maximumCutPasses_; 

2697  currentPassNumber_ = rhs.currentPassNumber_; 

2698  intParam_[CbcMaxNumNode] = rhs.intParam_[CbcMaxNumNode]; 

2699  intParam_[CbcMaxNumSol] = rhs.intParam_[CbcMaxNumSol]; 

2700  intParam_[CbcFathomDiscipline] = rhs.intParam_[CbcFathomDiscipline]; 

2701  dblParam_[CbcIntegerTolerance] = rhs.dblParam_[CbcIntegerTolerance]; 

2702  dblParam_[CbcInfeasibilityWeight] = rhs.dblParam_[CbcInfeasibilityWeight]; 

2703  dblParam_[CbcCutoffIncrement] = rhs.dblParam_[CbcCutoffIncrement]; 

2704  dblParam_[CbcAllowableGap] = rhs.dblParam_[CbcAllowableGap]; 

2705  dblParam_[CbcAllowableFractionGap] = rhs.dblParam_[CbcAllowableFractionGap]; 

2706  dblParam_[CbcMaximumSeconds] = rhs.dblParam_[CbcMaximumSeconds]; 

2707  dblParam_[CbcCurrentCutoff] = rhs.dblParam_[CbcCurrentCutoff]; 

2708  dblParam_[CbcOptimizationDirection] = rhs.dblParam_[CbcOptimizationDirection]; 

2709  dblParam_[CbcCurrentObjectiveValue] = rhs.dblParam_[CbcCurrentObjectiveValue]; 

2710  dblParam_[CbcCurrentMinimizationObjectiveValue] = rhs.dblParam_[CbcCurrentMinimizationObjectiveValue]; 

2711  dblParam_[CbcStartSeconds] = dblParam_[CbcStartSeconds]; // will be overwritten hopefully 

2712  globalCuts_ = rhs.globalCuts_; 

2713  int i; 

2714  for (i=0;i<numberCutGenerators_;i++) { 

2715  delete generator_[i]; 

2716  delete virginGenerator_[i]; 

2717  } 

2718  delete [] generator_; 

2719  delete [] virginGenerator_; 

2720  delete [] heuristic_; 

2721  maximumWhich_ = rhs.maximumWhich_; 

2722  delete [] whichGenerator_; 

2723  whichGenerator_ = NULL; 

2724  if (maximumWhich_&&rhs.whichGenerator_) 

2725  whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_); 

2726  for (i=0;i<maximumStatistics_;i++) 

2727  delete statistics_[i]; 

2728  delete [] statistics_; 

2729  maximumStatistics_ = 0; 

2730  statistics_ = NULL; 

2731  numberFixedAtRoot_ = rhs.numberFixedAtRoot_; 

2732  numberFixedNow_ = rhs.numberFixedNow_; 

2733  stoppedOnGap_ = rhs.stoppedOnGap_; 

2734  eventHappened_ = rhs.eventHappened_; 

2735  numberLongStrong_ = rhs.numberLongStrong_; 

2736  numberOldActiveCuts_ = rhs.numberOldActiveCuts_; 

2737  numberNewCuts_ = rhs.numberNewCuts_; 

2738  resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_; 

2739  sizeMiniTree_ = rhs.sizeMiniTree_; 

2740  searchStrategy_ = rhs.searchStrategy_; 

2741  numberStrongIterations_ = rhs.numberStrongIterations_; 

2742  strongInfo_[0]=rhs.strongInfo_[0]; 

2743  strongInfo_[1]=rhs.strongInfo_[1]; 

2744  strongInfo_[2]=rhs.strongInfo_[2]; 

2745  solverCharacteristics_ = NULL; 

2746  lastHeuristic_ = NULL; 

2747  numberCutGenerators_ = rhs.numberCutGenerators_; 

2748  if (numberCutGenerators_) { 

2749  generator_ = new CbcCutGenerator * [numberCutGenerators_]; 

2750  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_]; 

2751  int i; 

2752  for (i=0;i<numberCutGenerators_;i++) { 

2753  generator_[i]=new CbcCutGenerator(*rhs.generator_[i]); 

2754  virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]); 

2755  } 

2756  } else { 

2757  generator_=NULL; 

2758  virginGenerator_=NULL; 

2759  } 

2760  numberHeuristics_ = rhs.numberHeuristics_; 

2761  if (numberHeuristics_) { 

2762  heuristic_ = new CbcHeuristic * [numberHeuristics_]; 

2763  memcpy(heuristic_,rhs.heuristic_, 

2764  numberHeuristics_*sizeof(CbcHeuristic *)); 

2765  } else { 

2766  heuristic_=NULL; 

2767  } 

2768  lastHeuristic_ = NULL; 

2769  if (eventHandler_) 

2770  delete eventHandler_ ; 

2771  if (rhs.eventHandler_) 

2772  { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; } 

2773  else 

2774  { eventHandler_ = NULL ; } 

2775  for (i=0;i<numberObjects_;i++) 

2776  delete object_[i]; 

2777  delete [] object_; 

2778  numberObjects_=rhs.numberObjects_; 

2779  if (numberObjects_) { 

2780  object_ = new CbcObject * [numberObjects_]; 

2781  int i; 

2782  for (i=0;i<numberObjects_;i++) 

2783  object_[i]=(rhs.object_[i])>clone(); 

2784  } else { 

2785  object_=NULL; 

2786  } 

2787  delete [] originalColumns_; 

2788  if (rhs.originalColumns_) { 

2789  int numberColumns = rhs.getNumCols(); 

2790  originalColumns_= new int [numberColumns]; 

2791  memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int)); 

2792  } else { 

2793  originalColumns_=NULL; 

2794  } 

2795  nodeCompare_=rhs.nodeCompare_>clone(); 

2796  problemFeasibility_=rhs.problemFeasibility_>clone(); 

2797  delete tree_; 

2798  tree_= rhs.tree_>clone(); 

2799  branchingMethod_=rhs.branchingMethod_; 

2800  delete strategy_; 

2801  if (rhs.strategy_) 

2802  strategy_=rhs.strategy_>clone(); 

2803  else 

2804  strategy_=NULL; 

2805  parentModel_=rhs.parentModel_; 

2806  appData_=rhs.appData_; 

2807  

2808  delete [] integerVariable_; 

2809  numberIntegers_=rhs.numberIntegers_; 

2810  if (numberIntegers_) { 

2811  integerVariable_ = new int [numberIntegers_]; 

2812  memcpy(integerVariable_,rhs.integerVariable_, 

2813  numberIntegers_*sizeof(int)); 

2814  integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_>getNumCols()); 

2815  } else { 

2816  integerVariable_ = NULL; 

2817  integerInfo_=NULL; 

2818  } 

2819  if (rhs.hotstartSolution_) { 

2820  int numberColumns = solver_>getNumCols(); 

2821  hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns); 

2822  hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns); 

2823  } else { 

2824  hotstartSolution_ = NULL; 

2825  hotstartPriorities_ =NULL; 

2826  } 

2827  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_; 

2828  maximumNumberCuts_=rhs.maximumNumberCuts_; 

2829  phase_ = rhs.phase_; 

2830  currentNumberCuts_=rhs.currentNumberCuts_; 

2831  maximumDepth_= rhs.maximumDepth_; 

2832  delete [] addedCuts_; 

2833  delete [] walkback_; 

2834  // These are only used as temporary arrays so need not be filled 

2835  if (maximumNumberCuts_) { 

2836  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

2837  } else { 

2838  addedCuts_ = NULL; 

2839  } 

2840  nextRowCut_ = NULL; 

2841  currentNode_ = NULL; 

2842  if (maximumDepth_) 

2843  walkback_ = new CbcNodeInfo * [maximumDepth_]; 

2844  else 

2845  walkback_ = NULL; 

2846  synchronizeModel(); 

2847  cbcColLower_ = NULL; 

2848  cbcColUpper_ = NULL; 

2849  cbcRowLower_ = NULL; 

2850  cbcRowUpper_ = NULL; 

2851  cbcColSolution_ = NULL; 

2852  cbcRowPrice_ = NULL; 

2853  cbcReducedCost_ = NULL; 

2854  cbcRowActivity_ = NULL; 

2855  } 

2856  return *this; 

2857  } 

2858  

2859  // Destructor 

2860  CbcModel::~CbcModel () 

2861  { 

2862  if (defaultHandler_) { 

2863  delete handler_; 

2864  handler_ = NULL; 

2865  } 

2866  delete tree_; 

2867  tree_=NULL; 

2868  if (ourSolver_) delete solver_; 

2869  gutsOfDestructor(); 

2870  delete eventHandler_ ; 

2871  eventHandler_ = NULL ; 

2872  } 

2873  // Clears out as much as possible (except solver) 

2874  void 

2875  CbcModel::gutsOfDestructor() 

2876  { 

2877  delete referenceSolver_; 

2878  referenceSolver_=NULL; 

2879  int i; 

2880  for (i=0;i<numberCutGenerators_;i++) { 

2881  delete generator_[i]; 

2882  delete virginGenerator_[i]; 

2883  } 

2884  delete [] generator_; 

2885  delete [] virginGenerator_; 

2886  generator_=NULL; 

2887  virginGenerator_=NULL; 

2888  for (i=0;i<numberHeuristics_;i++) 

2889  delete heuristic_[i]; 

2890  delete [] heuristic_; 

2891  heuristic_=NULL; 

2892  delete nodeCompare_; 

2893  nodeCompare_=NULL; 

2894  delete problemFeasibility_; 

2895  problemFeasibility_=NULL; 

2896  delete [] originalColumns_; 

2897  originalColumns_=NULL; 

2898  delete strategy_; 

2899  gutsOfDestructor2(); 

2900  } 

2901  // Clears out enough to reset CbcModel 

2902  void 

2903  CbcModel::gutsOfDestructor2() 

2904  { 

2905  delete [] integerInfo_; 

2906  integerInfo_=NULL; 

2907  delete [] integerVariable_; 

2908  integerVariable_=NULL; 

2909  int i; 

2910  for (i=0;i<numberObjects_;i++) 

2911  delete object_[i]; 

2912  delete [] object_; 

2913  object_=NULL; 

2914  numberIntegers_=0; 

2915  numberObjects_=0; 

2916  delete emptyWarmStart_ ; 

2917  emptyWarmStart_ =NULL; 

2918  delete continuousSolver_; 

2919  continuousSolver_=NULL; 

2920  delete [] bestSolution_; 

2921  bestSolution_=NULL; 

2922  delete [] currentSolution_; 

2923  currentSolution_=NULL; 

2924  delete [] continuousSolution_; 

2925  continuousSolution_=NULL; 

2926  solverCharacteristics_=NULL; 

2927  delete [] usedInSolution_; 

2928  usedInSolution_ = NULL; 

2929  testSolution_=NULL; 

2930  lastHeuristic_ = NULL; 

2931  delete [] addedCuts_; 

2932  addedCuts_=NULL; 

2933  nextRowCut_ = NULL; 

2934  currentNode_ = NULL; 

2935  delete [] walkback_; 

2936  walkback_=NULL; 

2937  delete [] whichGenerator_; 

2938  whichGenerator_ = NULL; 

2939  for (i=0;i<maximumStatistics_;i++) 

2940  delete statistics_[i]; 

2941  delete [] statistics_; 

2942  statistics_=NULL; 

2943  maximumStatistics_=0; 

2944  delete [] analyzeResults_; 

2945  analyzeResults_=NULL; 

2946  // Below here is whatever consensus is 

2947  ourSolver_=true; 

2948  bestObjective_=COIN_DBL_MAX; 

2949  bestPossibleObjective_=COIN_DBL_MAX; 

2950  sumChangeObjective1_=0.0; 

2951  sumChangeObjective2_=0.0; 

2952  numberSolutions_=0; 

2953  stateOfSearch_=0; 

2954  delete [] hotstartSolution_; 

2955  hotstartSolution_=NULL; 

2956  delete [] hotstartPriorities_; 

2957  hotstartPriorities_=NULL; 

2958  numberHeuristicSolutions_=0; 

2959  numberNodes_=0; 

2960  numberNodes2_=0; 

2961  numberIterations_=0; 

2962  status_=1; 

2963  secondaryStatus_=1; 

2964  maximumNumberCuts_=0; 

2965  phase_=0; 

2966  currentNumberCuts_=0; 

2967  maximumDepth_=0; 

2968  nextRowCut_=NULL; 

2969  currentNode_=NULL; 

2970  // clear out tree 

2971  if (tree_&&tree_>size()) 

2972  tree_>cleanTree(this, 1.0e100,bestPossibleObjective_) ; 

2973  subTreeModel_=NULL; 

2974  numberStoppedSubTrees_=0; 

2975  numberInfeasibleNodes_=0; 

2976  numberGlobalViolations_=0; 

2977  continuousObjective_=0.0; 

2978  originalContinuousObjective_=0.0; 

2979  continuousInfeasibilities_=0; 

2980  numberFixedAtRoot_=0; 

2981  numberFixedNow_=0; 

2982  stoppedOnGap_=false; 

2983  eventHappened_=false; 

2984  numberLongStrong_=0; 

2985  numberOldActiveCuts_=0; 

2986  numberNewCuts_=0; 

2987  searchStrategy_=1; 

2988  numberStrongIterations_=0; 

2989  // Parameters which need to be reset 

2990  dblParam_[CbcCutoffIncrement] = 1e5; 

2991  dblParam_[CbcCurrentCutoff] = 1.0e100; 

2992  dblParam_[CbcCurrentObjectiveValue] = 1.0e100; 

2993  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100; 

2994  } 

2995  // Save a copy of the current solver so can be reset to 

2996  void 

2997  CbcModel::saveReferenceSolver() 

2998  { 

2999  delete referenceSolver_; 

3000  referenceSolver_= solver_>clone(); 

3001  } 

3002  

3003  // Uses a copy of reference solver to be current solver 

3004  void 

3005  CbcModel::resetToReferenceSolver() 

3006  { 

3007  delete solver_; 

3008  solver_ = referenceSolver_>clone(); 

3009  // clear many things 

3010  gutsOfDestructor2(); 

3011  // Reset cutoff 

3012  // Solvers know about direction 

3013  double direction = solver_>getObjSense(); 

3014  double value; 

3015  solver_>getDblParam(OsiDualObjectiveLimit,value); 

3016  setCutoff(value*direction); 

3017  } 

3018  

3019  // Are there a numerical difficulties? 

3020  bool 

3021  CbcModel::isAbandoned() const 

3022  { 

3023  return status_ == 2; 

3024  } 

3025  // Is optimality proven? 

3026  bool 

3027  CbcModel::isProvenOptimal() const 

3028  { 

3029  if (!status_ && bestObjective_<1.0e30) 

3030  return true; 

3031  else 

3032  return false; 

3033  } 

3034  // Is infeasiblity proven (or none better than cutoff)? 

3035  bool 

3036  CbcModel::isProvenInfeasible() const 

3037  { 

3038  if (!status_ && bestObjective_>=1.0e30) 

3039  return true; 

3040  else 

3041  return false; 

3042  } 

3043  // Node limit reached? 

3044  bool 

3045  CbcModel::isNodeLimitReached() const 

3046  { 

3047  return numberNodes_ >= intParam_[CbcMaxNumNode]; 

3048  } 

3049  // Time limit reached? 

3050  bool 

3051  CbcModel::isSecondsLimitReached() const 

3052  { 

3053  if (status_==1&&secondaryStatus_==4) 

3054  return true; 

3055  else 

3056  return false; 

3057  } 

3058  // Solution limit reached? 

3059  bool 

3060  CbcModel::isSolutionLimitReached() const 

3061  { 

3062  return numberSolutions_ >= intParam_[CbcMaxNumSol]; 

3063  } 

3064  // Set language 

3065  void 

3066  CbcModel::newLanguage(CoinMessages::Language language) 

3067  { 

3068  messages_ = CbcMessage(language); 

3069  } 

3070  void 

3071  CbcModel::setNumberStrong(int number) 

3072  { 

3073  if (number<0) 

3074  numberStrong_=0; 

3075  else 

3076  numberStrong_=number; 

3077  } 

3078  void 

3079  CbcModel::setNumberBeforeTrust(int number) 

3080  { 

3081  if (number<1) { 

3082  numberBeforeTrust_=0; 

3083  } else { 

3084  numberBeforeTrust_=number; 

3085  //numberStrong_ = CoinMax(numberStrong_,1); 

3086  } 

3087  } 

3088  void 

3089  CbcModel::setNumberPenalties(int number) 

3090  { 

3091  if (number<=0) { 

3092  numberPenalties_=0; 

3093  } else { 

3094  numberPenalties_=number; 

3095  } 

3096  } 

3097  void 

3098  CbcModel::setPenaltyScaleFactor(double value) 

3099  { 

3100  if (value<=0) { 

3101  penaltyScaleFactor_=3.0; 

3102  } else { 

3103  penaltyScaleFactor_=value; 

3104  } 

3105  } 

3106  void 

3107  CbcModel::setHowOftenGlobalScan(int number) 

3108  { 

3109  if (number<1) 

3110  howOftenGlobalScan_=0; 

3111  else 

3112  howOftenGlobalScan_=number; 

3113  } 

3114  

3115  // Add one generator 

3116  void 

3117  CbcModel::addCutGenerator(CglCutGenerator * generator, 

3118  int howOften, const char * name, 

3119  bool normal, bool atSolution, 

3120  bool whenInfeasible,int howOftenInSub, 

3121  int whatDepth, int whatDepthInSub) 

3122  { 

3123  CbcCutGenerator ** temp = generator_; 

3124  generator_ = new CbcCutGenerator * [numberCutGenerators_+1]; 

3125  memcpy(generator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *)); 

3126  delete[] temp ; 

3127  generator_[numberCutGenerators_]= 

3128  new CbcCutGenerator(this,generator, howOften, name, 

3129  normal,atSolution,whenInfeasible,howOftenInSub, 

3130  whatDepth, whatDepthInSub); 

3131  // and before any cahnges 

3132  temp = virginGenerator_; 

3133  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_+1]; 

3134  memcpy(virginGenerator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *)); 

3135  delete[] temp ; 

3136  virginGenerator_[numberCutGenerators_++]= 

3137  new CbcCutGenerator(this,generator, howOften, name, 

3138  normal,atSolution,whenInfeasible,howOftenInSub, 

3139  whatDepth, whatDepthInSub); 

3140  

3141  } 

3142  // Add one heuristic 

3143  void 

3144  CbcModel::addHeuristic(CbcHeuristic * generator) 

3145  { 

3146  CbcHeuristic ** temp = heuristic_; 

3147  heuristic_ = new CbcHeuristic * [numberHeuristics_+1]; 

3148  memcpy(heuristic_,temp,numberHeuristics_*sizeof(CbcHeuristic *)); 

3149  delete [] temp; 

3150  heuristic_[numberHeuristics_++]=generator>clone(); 

3151  } 

3152  

3153  /* 

3154  The last subproblem handled by the solver is not necessarily related to the 

3155  one being recreated, so the first action is to remove all cuts from the 

3156  constraint system. Next, traverse the tree from node to the root to 

3157  determine the basis size required for this subproblem and create an empty 

3158  basis with the right capacity. Finally, traverse the tree from root to 

3159  node, adjusting bounds in the constraint system, adjusting the basis, and 

3160  collecting the cuts that must be added to the constraint system. 

3161  applyToModel does the heavy lifting. 

3162  

3163  addCuts1 is used in contexts where all that's desired is the list of cuts: 

3164  the node is already fathomed, and we're collecting cuts so that we can 

3165  adjust reference counts as we prune nodes. Arguably the two functions 

3166  should be separated. The culprit is applyToModel, which performs cut 

3167  collection and model adjustment. 

3168  

3169  Certainly in the contexts where all we need is a list of cuts, there's no 

3170  point in passing in a valid basis  an empty basis will do just fine. 

3171  */ 

3172  void CbcModel::addCuts1 (CbcNode * node, CoinWarmStartBasis *&lastws) 

3173  { int i; 

3174  int nNode=0; 

3175  int numberColumns = getNumCols(); 

3176  CbcNodeInfo * nodeInfo = node>nodeInfo(); 

3177  

3178  /* 

3179  Remove all cuts from the constraint system. 

3180  (original comment includes ``see note below for later efficiency'', but 

3181  the reference isn't clear to me). 

3182  */ 

3183  int currentNumberCuts = solver_>getNumRows()numberRowsAtContinuous_; 

3184  int *which = new int[currentNumberCuts]; 

3185  for (i = 0 ; i < currentNumberCuts ; i++) 

3186  which[i] = i+numberRowsAtContinuous_; 

3187  solver_>deleteRows(currentNumberCuts,which); 

3188  delete [] which; 

3189  /* 

3190  Accumulate the path from node to the root in walkback_, and accumulate a 

3191  cut count in currentNumberCuts. 

3192  

3193  original comment: when working then just unwind until where new node joins 

3194  old node (for cuts?) 

3195  */ 

3196  currentNumberCuts = 0; 

3197  while (nodeInfo) { 

3198  //printf("nNode = %d, nodeInfo = %x\n",nNode,nodeInfo); 

3199  walkback_[nNode++]=nodeInfo; 

3200  currentNumberCuts += nodeInfo>numberCuts() ; 

3201  nodeInfo = nodeInfo>parent() ; 

3202  if (nNode==maximumDepth_) { 

3203  maximumDepth_ *= 2; 

3204  CbcNodeInfo ** temp = new CbcNodeInfo * [maximumDepth_]; 

3205  for (i=0;i<nNode;i++) 

3206  temp[i] = walkback_[i]; 

3207  delete [] walkback_; 

3208  walkback_ = temp; 

3209  } 

3210  } 

3211  /* 

3212  Create an empty basis with sufficient capacity for the constraint system 

3213  we'll construct: original system plus cuts. Make sure we have capacity to 

3214  record those cuts in addedCuts_. 

3215  

3216  The method of adjusting the basis at a FullNodeInfo object (the root, for 

3217  example) is to use a copy constructor to duplicate the basis held in the 

3218  nodeInfo, then resize it and return the new basis object. Guaranteed, 

3219  lastws will point to a different basis when it returns. We pass in a basis 

3220  because we need the parameter to return the allocated basis, and it's an 

3221  easy way to pass in the size. But we take a hit for memory allocation. 

3222  */ 

3223  currentNumberCuts_=currentNumberCuts; 

3224  if (currentNumberCuts >= maximumNumberCuts_) { 

3225  maximumNumberCuts_ = currentNumberCuts; 

3226  delete [] addedCuts_; 

3227  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

3228  } 

3229  lastws>setSize(numberColumns,numberRowsAtContinuous_+currentNumberCuts); 

3230  /* 

3231  This last bit of code traverses the path collected in walkback_ from the 

3232  root back to node. At the end of the loop, 

3233  * lastws will be an appropriate basis for node; 

3234  * variable bounds in the constraint system will be set to be correct for 

3235  node; and 

3236  * addedCuts_ will be set to a list of cuts that need to be added to the 

3237  constraint system at node. 

3238  applyToModel does all the heavy lifting. 

3239  */ 

3240  currentNumberCuts=0; 

3241  while (nNode) { 

3242  nNode; 

3243  walkback_[nNode]>applyToModel(this,lastws,addedCuts_,currentNumberCuts); 

3244  } 

3245  } 

3246  

3247  /* 

3248  adjustCuts might be a better name: If the node is feasible, we sift through 

3249  the cuts collected by addCuts1, add the ones that are tight and omit the 

3250  ones that are loose. If the node is infeasible, we just adjust the 

3251  reference counts to reflect that we're about to prune this node and its 

3252  descendants. 

3253  */ 

3254  int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws,bool canFix) 

3255  { 

3256  /* 

3257  addCuts1 performs step 1 of restoring the subproblem at this node; see the 

3258  comments there. 

3259  */ 

3260  addCuts1(node,lastws); 

3261  int i; 

3262  int numberColumns = getNumCols(); 

3263  CbcNodeInfo * nodeInfo = node>nodeInfo(); 

3264  double cutoff = getCutoff() ; 

3265  int currentNumberCuts=currentNumberCuts_; 

3266  if (canFix) { 

3267  bool feasible=true; 

3268  const double *lower = solver_>getColLower() ; 

3269  const double *upper = solver_>getColUpper() ; 

3270  double * newLower = analyzeResults_; 

3271  double * objLower = newLower+numberIntegers_; 

3272  double * newUpper = objLower+numberIntegers_; 

3273  double * objUpper = newUpper+numberIntegers_; 

3274  int n=0; 

3275  for (i=0;i<numberIntegers_;i++) { 

3276  int iColumn = integerVariable_[i]; 

3277  bool changed=false; 

3278  double lo = 0.0; 

3279  double up = 0.0; 

3280  if (objLower[i]>cutoff) { 

3281  lo = lower[iColumn]; 

3282  up = upper[iColumn]; 

3283  if (lo<newLower[i]) { 

3284  lo = newLower[i]; 

3285  solver_>setColLower(iColumn,lo); 

3286  changed=true; 

3287  n++; 

3288  } 

3289  if (objUpper[i]>cutoff) { 

3290  if (up>newUpper[i]) { 

3291  up = newUpper[i]; 

3292  solver_>setColUpper(iColumn,up); 

3293  changed=true; 

3294  n++; 

3295  } 

3296  } 

3297  } else if (objUpper[i]>cutoff) { 

3298  lo = lower[iColumn]; 

3299  up = upper[iColumn]; 

3300  if (up>newUpper[i]) { 

3301  up = newUpper[i]; 

3302  solver_>setColUpper(iColumn,up); 

3303  changed=true; 

3304  n++; 

3305  } 

3306  } 

3307  if (changed&&lo>up) { 

3308  feasible=false; 

3309  break; 

3310  } 

3311  } 

3312  if (!feasible) { 

3313  printf("analysis says node infeas\n"); 

3314  cutoff=COIN_DBL_MAX; 

3315  } 

3316  } 

3317  /* 

3318  If the node can't be fathomed by bound, reinstall tight cuts in the 

3319  constraint system. Even if there are no cuts, we'll want to set the 

3320  reconstructed basis in the solver. 

3321  */ 

3322  if (node>objectiveValue() < cutoff) 

3323  { 

3324  # ifdef CBC_CHECK_BASIS 

3325  printf("addCuts: expanded basis; rows %d+%d\n", 

3326  numberRowsAtContinuous_,currentNumberCuts); 

3327  lastws>print(); 

3328  # endif 

3329  /* 

3330  Adjust the basis and constraint system so that we retain only active cuts. 

3331  There are three steps: 

3332  1) Scan the basis. Sort the cuts into effective cuts to be kept and 

3333  loose cuts to be dropped. 

3334  2) Drop the loose cuts and resize the basis to fit. 

3335  3) Install the tight cuts in the constraint system (applyRowCuts) and 

3336  and install the basis (setWarmStart). 

3337  Use of compressRows conveys we're compressing the basis and not just 

3338  tweaking the artificialStatus_ array. 

3339  */ 

3340  if (currentNumberCuts > 0) { 

3341  int numberToAdd = 0; 

3342  const OsiRowCut **addCuts; 

3343  int numberToDrop = 0 ; 

3344  int *cutsToDrop ; 

3345  addCuts = new const OsiRowCut* [currentNumberCuts]; 

3346  cutsToDrop = new int[currentNumberCuts] ; 

3347  int nxrow = lastws>getNumArtificial(); 

3348  for (i=0;i<currentNumberCuts;i++) { 

3349  assert (i+numberRowsAtContinuous_<nxrow); 

3350  CoinWarmStartBasis::Status status = 

3351  lastws>getArtifStatus(i+numberRowsAtContinuous_); 

3352  if (addedCuts_[i] && 

3353  (status != CoinWarmStartBasis::basic  

3354  addedCuts_[i]>effectiveness()==COIN_DBL_MAX)) { 

3355  # ifdef CHECK_CUT_COUNTS 

3356  printf("Using cut %d %x as row %d\n",i,addedCuts_[i], 

3357  numberRowsAtContinuous_+numberToAdd); 

3358  # endif 

3359  addCuts[numberToAdd++] = new OsiRowCut(*addedCuts_[i]); 

3360  } else { 

3361  # ifdef CHECK_CUT_COUNTS 

3362  printf("Dropping cut %d %x\n",i,addedCuts_[i]); 

3363  # endif 

3364  addedCuts_[i]=NULL; 

3365  cutsToDrop[numberToDrop++] = numberRowsAtContinuous_+i ; 

3366  } 

3367  } 

3368  int numberRowsNow=numberRowsAtContinuous_+numberToAdd; 

3369  lastws>compressRows(numberToDrop,cutsToDrop) ; 

3370  lastws>resize(numberRowsNow,numberColumns); 

3371  solver_>applyRowCuts(numberToAdd,addCuts); 

3372  # ifdef CBC_CHECK_BASIS 

3373  printf("addCuts: stripped basis; rows %d + %d\n", 

3374  numberRowsAtContinuous_,numberToAdd); 

3375  lastws>print(); 

3376  # endif 

3377  for (i=0;i<numberToAdd;i++) 

3378  delete addCuts[i]; 

3379  delete [] addCuts; 

3380  delete [] cutsToDrop ; 

3381  } 

3382  /* 

3383  Set the basis in the solver. 

3384  */ 

3385  solver_>setWarmStart(lastws); 

3386  #if 0 

3387  if ((numberNodes_%printFrequency_)==0) { 

3388  printf("Objective %g, depth %d, unsatisfied %d\n", 

3389  node>objectiveValue(), 

3390  node>depth(),node>numberUnsatisfied()); 

3391  } 

3392  #endif 

3393  /* 

3394  Clean up and we're out of here. 

3395  */ 

3396  numberNodes_++; 

3397  return 0; 

3398  } 

3399  /* 

3400  This node has been fathomed by bound as we try to revive it out of the live 

3401  set. Adjust the cut reference counts to reflect that we no longer need to 

3402  explore the remaining branch arms, hence they will no longer reference any 

3403  cuts. Cuts whose reference count falls to zero are deleted. 

3404  */ 

3405  else 

3406  { int i; 

3407  int numberLeft = nodeInfo>numberBranchesLeft(); 

3408  for (i = 0 ; i < currentNumberCuts ; i++) 

3409  { if (addedCuts_[i]) 

3410  { if (!addedCuts_[i]>decrement(numberLeft)) 

3411  { delete addedCuts_[i]; 

3412  addedCuts_[i] = NULL; } } } 

3413  return 1 ; } 

3414  } 

3415  

3416  

3417  /* 

3418  Perform reduced cost fixing on integer variables. 

3419  

3420  The variables in question are already nonbasic at bound. We're just nailing 

3421  down the current situation. 

3422  */ 

3423  

3424  int CbcModel::reducedCostFix () 

3425  

3426  { 

3427  if(!solverCharacteristics_>reducedCostsAccurate()) 

3428  return 0; //NLP 

3429  double cutoff = getCutoff() ; 

3430  double direction = solver_>getObjSense() ; 

3431  double gap = cutoff  solver_>getObjValue()*direction ; 

3432  double tolerance; 

3433  solver_>getDblParam(OsiDualTolerance,tolerance) ; 

3434  if (gap<=0.0) 

3435  return 0; 

3436  gap += 100.0*tolerance; 

3437  double integerTolerance = getDblParam(CbcIntegerTolerance) ; 

3438  

3439  const double *lower = solver_>getColLower() ; 

3440  const double *upper = solver_>getColUpper() ; 

3441  const double *solution = solver_>getColSolution() ; 

3442  const double *reducedCost = solver_>getReducedCost() ; 

3443  

3444  int numberFixed = 0 ; 

3445  

3446  # ifdef CBC_USE_CLP 

3447  OsiClpSolverInterface * clpSolver 

3448  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

3449  ClpSimplex * clpSimplex=NULL; 

3450  if (clpSolver) 

3451  clpSimplex = clpSolver>getModelPtr(); 

3452  # endif 

3453  for (int i = 0 ; i < numberIntegers_ ; i++) 

3454  { int iColumn = integerVariable_[i] ; 

3455  double djValue = direction*reducedCost[iColumn] ; 

3456  if (upper[iColumn]lower[iColumn] > integerTolerance) 

3457  { if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap) 

3458  { solver_>setColUpper(iColumn,lower[iColumn]) ; 

3459  #ifdef CBC_USE_CLP 

3460  if (clpSimplex) 

3461  assert (clpSimplex>getColumnStatus(iColumn)==ClpSimplex::atLowerBound); 

3462  #endif 

3463  numberFixed++ ; } 

3464  else 

3465  if (solution[iColumn] > upper[iColumn]integerTolerance && djValue > gap) 

3466  { solver_>setColLower(iColumn,upper[iColumn]) ; 

3467  #ifdef CBC_USE_CLP 

3468  if (clpSimplex) 

3469  assert (clpSimplex>getColumnStatus(iColumn)==ClpSimplex::atUpperBound); 

3470  #endif 

3471  numberFixed++ ; } } } 

3472  

3473  return numberFixed; } 

3474  

3475  // Collect coding to replace whichGenerator 

3476  void 

3477  CbcModel::resizeWhichGenerator(int numberNow, int numberAfter) 

3478  { 

3479  if (numberAfter > maximumWhich_) { 

3480  maximumWhich_ = CoinMax(maximumWhich_*2+100,numberAfter) ; 

3481  int * temp = new int[2*maximumWhich_] ; 

3482  memcpy(temp,whichGenerator_,numberNow*sizeof(int)) ; 

3483  delete [] whichGenerator_ ; 

3484  whichGenerator_ = temp ; 

3485  memset(whichGenerator_+numberNow,0,(maximumWhich_numberNow)*sizeof(int)); 

3486  } 

3487  } 

3488  

3489  /** Solve the model using cuts 

3490  

3491  This version takes off redundant cuts from node. 

3492  Returns true if feasible. 

3493  

3494  \todo 

3495  Why do I need to resolve the problem? What has been done between the last 

3496  relaxation and calling solveWithCuts? 

3497  

3498  If numberTries == 0 then user did not want any cuts. 

3499  */ 

3500  

3501  bool 

3502  CbcModel::solveWithCuts (OsiCuts &cuts, int numberTries, CbcNode *node) 

3503  /* 

3504  Parameters: 

3505  numberTries: (i) the maximum number of iterations for this round of cut 

3506  generation; if negative then we don't mind if drop is tiny. 

3507  

3508  cuts: (o) all cuts generated in this round of cut generation 

3509  

3510  node: (i) So we can update dynamic pseudo costs 

3511  */ 

3512  

3513  

3514  { 

3515  bool feasible = true ; 

3516  int lastNumberCuts = 0 ; 

3517  double lastObjective = 1.0e100 ; 

3518  int violated = 0 ; 

3519  int numberRowsAtStart = solver_>getNumRows() ; 

3520  //printf("solver had %d rows\n",numberRowsAtStart); 

3521  int numberColumns = solver_>getNumCols() ; 

3522  CoinBigIndex numberElementsAtStart = solver_>getNumElements(); 

3523  

3524  numberOldActiveCuts_ = numberRowsAtStartnumberRowsAtContinuous_ ; 

3525  numberNewCuts_ = 0 ; 

3526  

3527  bool onOptimalPath = false ; 

3528  const OsiRowCutDebugger *debugger = NULL; 

3529  if ((specialOptions_&1)!=0) { 

3530  /* 

3531  See OsiRowCutDebugger for details. In a nutshell, make sure that current 

3532  variable values do not conflict with a known optimal solution. (Obviously 

3533  this can be fooled when there are multiple solutions.) 

3534  */ 

3535  debugger = solver_>getRowCutDebugger() ; 

3536  if (debugger) 

3537  onOptimalPath = (debugger>onOptimalPath(*solver_)) ; 

3538  } 

3539  OsiCuts slackCuts; 

3540  /* 

3541  Resolve the problem. If we've lost feasibility, might as well bail out right 

3542  after the debug stuff. The resolve will also refresh cached copies of the 

3543  solver solution (cbcColLower_, ...) held by CbcModel. 

3544  */ 

3545  double objectiveValue = solver_>getObjValue()*solver_>getObjSense(); 

3546  if (node) 

3547  objectiveValue= node>objectiveValue(); 

3548  int returnCode = resolve(node ? node>nodeInfo() : NULL,1); 

3549  if (node&&!node>nodeInfo()>numberBranchesLeft()) 

3550  node>nodeInfo()>allBranchesGone(); // can clean up 

3551  feasible = returnCode != 0 ; 

3552  if (returnCode<0) 

3553  numberTries=0; 

3554  if (problemFeasibility_>feasible(this,0)<0) { 

3555  feasible=false; // pretend infeasible 

3556  } 

3557  

3558  // Update branching information if wanted 

3559  if(node &&branchingMethod_) 

3560  branchingMethod_>updateInformation(solver_,node); 

3561  

3562  #ifdef CBC_DEBUG 

3563  if (feasible) 

3564  { printf("Obj value %g (%s) %d rows\n",solver_>getObjValue(), 

3565  (solver_>isProvenOptimal())?"proven":"unproven", 

3566  solver_>getNumRows()) ; } 

3567  

3568  else 

3569  { printf("Infeasible %d rows\n",solver_>getNumRows()) ; } 

3570  #endif 

3571  if ((specialOptions_&1)!=0) { 

3572  /* 

3573  If the RowCutDebugger said we were compatible with the optimal solution, 

3574  and now we're suddenly infeasible, we might be confused. Then again, we 

3575  may have fathomed by bound, heading for a rediscovery of an optimal solution. 

3576  */ 

3577  if (onOptimalPath && !solver_>isDualObjectiveLimitReached()) { 

3578  if (!feasible) { 

3579  solver_>writeMps("infeas"); 

3580  CoinWarmStartBasis *slack = 

3581  dynamic_cast<CoinWarmStartBasis *>(solver_>getEmptyWarmStart()) ; 

3582  solver_>setWarmStart(slack); 

3583  delete slack ; 

3584  solver_>setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ; 

3585  solver_>initialSolve(); 

3586  } 

3587  assert(feasible) ; 

3588  } 

3589  } 

3590  

3591  if (!feasible) { 

3592  numberInfeasibleNodes_++; 

3593  return (false) ; 

3594  } 

3595  sumChangeObjective1_ += solver_>getObjValue()*solver_>getObjSense() 

3596   objectiveValue ; 

3597  if ( CoinCpuTime()dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] ) 

3598  numberTries=0; // exit 

3599  //if ((numberNodes_%100)==0) 

3600  //printf("XXa sum obj changed by %g\n",sumChangeObjective1_); 

3601  objectiveValue = solver_>getObjValue()*solver_>getObjSense(); 

3602  // Return at once if numberTries zero 

3603  if (!numberTries) { 

3604  cuts=OsiCuts(); 

3605  numberNewCuts_=0; 

3606  return true; 

3607  } 

3608  /* 

3609  Do reduced cost fixing. 

3610  */ 

3611  reducedCostFix() ; 

3612  /* 

3613  Set up for at most numberTries rounds of cut generation. If numberTries is 

3614  negative, we'll ignore the minimumDrop_ cutoff and keep generating cuts for 

3615  the specified number of rounds. 

3616  */ 

3617  double minimumDrop = minimumDrop_ ; 

3618  if (numberTries<0) 

3619  { numberTries = numberTries ; 

3620  minimumDrop = 1.0 ; } 

3621  /* 

3622  Is it time to scan the cuts in order to remove redundant cuts? If so, set 

3623  up to do it. 

3624  */ 

3625  # define SCANCUTS 100 

3626  int *countColumnCuts = NULL ; 

3627  // Always accumulate row cut counts 

3628  int * countRowCuts =new int[numberCutGenerators_+numberHeuristics_] ; 

3629  memset(countRowCuts,0, 

3630  (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ; 

3631  bool fullScan = false ; 

3632  if ((numberNodes_%SCANCUTS) == 0) 

3633  { fullScan = true ; 

3634  countColumnCuts = new int[numberCutGenerators_+numberHeuristics_] ; 

3635  memset(countColumnCuts,0, 

3636  (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ; } 

3637  

3638  double direction = solver_>getObjSense() ; 

3639  double startObjective = solver_>getObjValue()*direction ; 

3640  

3641  currentPassNumber_ = 0 ; 

3642  double primalTolerance = 1.0e7 ; 

3643  /* 

3644  Begin cut generation loop. Cuts generated during each iteration are 

3645  collected in theseCuts. The loop can be divided into four phases: 

3646  1) Prep: Fix variables using reduced cost. In the first iteration only, 

3647  consider scanning globalCuts_ and activating any applicable cuts. 

3648  2) Cut Generation: Call each generator and heuristic registered in the 

3649  generator_ and heuristic_ arrays. Newly generated global cuts are 

3650  copied to globalCuts_ at this time. 

3651  3) Cut Installation and Reoptimisation: Install column and row cuts in 

3652  the solver. Copy row cuts to cuts (parameter). Reoptimise. 

3653  4) Cut Purging: takeOffCuts() removes inactive cuts from the solver, and 

3654  does the necessary bookkeeping in the model. 

3655  */ 

3656  do 

3657  { currentPassNumber_++ ; 

3658  numberTries ; 

3659  OsiCuts theseCuts ; 

3660  /* 

3661  Scan previously generated global column and row cuts to see if any are 

3662  useful. 

3663  */ 

3664  int numberViolated=0; 

3665  if (currentPassNumber_ == 1 && howOftenGlobalScan_ > 0 && 

3666  (numberNodes_%howOftenGlobalScan_) == 0) 

3667  { int numberCuts = globalCuts_.sizeColCuts() ; 

3668  int i; 

3669  // possibly extend whichGenerator 

3670  resizeWhichGenerator(numberViolated, numberViolated+numberCuts); 

3671  for ( i = 0 ; i < numberCuts ; i++) 

3672  { const OsiColCut *thisCut = globalCuts_.colCutPtr(i) ; 

3673  if (thisCut>violated(cbcColSolution_)>primalTolerance) { 

3674  printf("Global cut added  violation %g\n", 

3675  thisCut>violated(cbcColSolution_)) ; 

3676  whichGenerator_[numberViolated++]=1; 

3677  theseCuts.insert(*thisCut) ; 

3678  } 

3679  } 

3680  numberCuts = globalCuts_.sizeRowCuts() ; 

3681  // possibly extend whichGenerator 

3682  resizeWhichGenerator(numberViolated, numberViolated+numberCuts); 

3683  for ( i = 0;i<numberCuts;i++) { 

3684  const OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ; 

3685  if (thisCut>violated(cbcColSolution_)>primalTolerance) { 

3686  //printf("Global cut added  violation %g\n", 

3687  // thisCut>violated(cbcColSolution_)) ; 

3688  whichGenerator_[numberViolated++]=1; 

3689  theseCuts.insert(*thisCut) ; 

3690  } 

3691  } 

3692  numberGlobalViolations_+=numberViolated; 

3693  } 

3694  /* 

3695  Generate new cuts (global and/or local) and/or apply heuristics. If 

3696  CglProbing is used, then it should be first as it can fix continuous 

3697  variables. 

3698  

3699  At present, CglProbing is the only case where generateCuts will return 

3700  true. generateCuts actually modifies variable bounds in the solver when 

3701  CglProbing indicates that it can fix a variable. Reoptimisation is required 

3702  to take full advantage. 

3703  

3704  The need to resolve here should only happen after a heuristic solution. 

3705  (Note default OSI implementation of optimalBasisIsAvailable always returns 

3706  false.) 

3707  */ 

3708  if (solverCharacteristics_>warmStart()&& 

3709  !solver_>optimalBasisIsAvailable()) { 

3710  //printf("XXXXYY no opt basis\n"); 

3711  resolve(node ? node>nodeInfo() : NULL,3); 

3712  } 

3713  if (nextRowCut_) { 

3714  // branch was a cut  add it 

3715  theseCuts.insert(*nextRowCut_); 

3716  if (handler_>logLevel()>1) 

3717  nextRowCut_>print(); 

3718  const OsiRowCut * cut=nextRowCut_; 

3719  double lb = cut>lb(); 

3720  double ub = cut>ub(); 

3721  int n=cut>row().getNumElements(); 

3722  const int * column = cut>row().getIndices(); 

3723  const double * element = cut>row().getElements(); 

3724  double sum=0.0; 

3725  for (int i=0;i<n;i++) { 

3726  int iColumn = column[i]; 

3727  double value = element[i]; 

3728  //if (cbcColSolution_[iColumn]>1.0e7) 

3729  //printf("value of %d is %g\n",iColumn,cbcColSolution_[iColumn]); 

3730  sum += value * cbcColSolution_[iColumn]; 

3731  } 

3732  delete nextRowCut_; 

3733  nextRowCut_=NULL; 

3734  if (handler_>logLevel()>1) 

3735  printf("applying branch cut, sum is %g, bounds %g %g\n",sum,lb,ub); 

3736  // possibly extend whichGenerator 

3737  resizeWhichGenerator(numberViolated, numberViolated+1); 

3738  // set whichgenerator (also serves as marker to say don't delete0 

3739  whichGenerator_[numberViolated++]=2; 

3740  } 

3741  double * newSolution = new double [numberColumns] ; 

3742  double heuristicValue = getCutoff() ; 

3743  int found = 1; // no solution found 

3744  

3745  for (int i = 0;i<numberCutGenerators_+numberHeuristics_;i++) { 

3746  int numberRowCutsBefore = theseCuts.sizeRowCuts() ; 

3747  int numberColumnCutsBefore = theseCuts.sizeColCuts() ; 

3748  if (i<numberCutGenerators_) { 

3749  bool generate = generator_[i]>normal(); 

3750  // skip if not optimal and should be (maybe a cut generator has fixed variables) 

3751  if (generator_[i]>needsOptimalBasis()&&!solver_>basisIsAvailable()) 

3752  generate=false; 

3753  if (generate) { 

3754  bool mustResolve = 

3755  generator_[i]>generateCuts(theseCuts,fullScan,node) ; 

3756  #ifdef CBC_DEBUG 

3757  { 

3758  int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 

3759  int k ; 

3760  for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 

3761  OsiRowCut thisCut = theseCuts.rowCut(k) ; 

3762  /* check size of elements. 

3763  We can allow smaller but this helps debug generators as it 

3764  is unsafe to have small elements */ 

3765  int n=thisCut.row().getNumElements(); 

3766  const int * column = thisCut.row().getIndices(); 

3767  const double * element = thisCut.row().getElements(); 

3768  //assert (n); 

3769  for (int i=0;i<n;i++) { 

3770  double value = element[i]; 

3771  assert(fabs(value)>1.0e12&&fabs(value)<1.0e20); 

3772  } 

3773  } 

3774  } 

3775  #endif 

3776  if (mustResolve) { 

3777  int returncode = resolve(node ? node>nodeInfo() : NULL,2); 

3778  feasible = returnCode != 0 ; 

3779  if (returncode<0) 

3780  numberTries=0; 

3781  if ((specialOptions_&1)!=0) { 

3782  debugger = solver_>getRowCutDebugger() ; 

3783  if (debugger) 

3784  onOptimalPath = (debugger>onOptimalPath(*solver_)) ; 

3785  else 

3786  onOptimalPath=false; 

3787  if (onOptimalPath && !solver_>isDualObjectiveLimitReached()) 

3788  assert(feasible) ; 

3789  } 

3790  if (!feasible) 

3791  break ; 

3792  } 

3793  } 

3794  } else { 

3795  // see if heuristic will do anything 

3796  double saveValue = heuristicValue ; 

3797  int ifSol = 

3798  heuristic_[inumberCutGenerators_]>solution(heuristicValue, 

3799  newSolution, 

3800  theseCuts) ; 

3801  if (ifSol>0) { 

3802  // better solution found 

3803  found = inumberCutGenerators_ ; 

3804  incrementUsed(newSolution); 

3805  } else if (ifSol<0) { 

3806  heuristicValue = saveValue ; 

3807  } 

3808  } 

3809  int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 

3810  int numberColumnCutsAfter = theseCuts.sizeColCuts() ; 

3811  

3812  if ((specialOptions_&1)!=0) { 

3813  if (onOptimalPath) { 

3814  int k ; 

3815  for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 

3816  OsiRowCut thisCut = theseCuts.rowCut(k) ; 

3817  if(debugger>invalidCut(thisCut)) { 

3818  solver_>writeMps("badCut"); 

3819  } 

3820  assert(!debugger>invalidCut(thisCut)) ; 

3821  } 

3822  } 

3823  } 

3824  /* 

3825  The cut generator/heuristic has done its thing, and maybe it generated some 

3826  cuts and/or a new solution. Do a bit of bookkeeping: load 

3827  whichGenerator[i] with the index of the generator responsible for a cut, 

3828  and place cuts flagged as global in the global cut pool for the model. 

3829  

3830  lastNumberCuts is the sum of cuts added in previous iterations; it's the 

3831  offset to the proper starting position in whichGenerator. 

3832  */ 

3833  int numberBefore = 

3834  numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ; 

3835  int numberAfter = 

3836  numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ; 

3837  // possibly extend whichGenerator 

3838  resizeWhichGenerator(numberBefore, numberAfter); 

3839  int j ; 

3840  if (fullScan) { 

3841  // counts 

3842  countColumnCuts[i] += numberColumnCutsAfternumberColumnCutsBefore ; 

3843  } 

3844  countRowCuts[i] += numberRowCutsAfternumberRowCutsBefore ; 

3845  

3846  for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) { 

3847  whichGenerator_[numberBefore++] = i ; 

3848  const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ; 

3849  if (thisCut>lb()>thisCut>ub()) 

3850  violated=2; // subproblem is infeasible 

3851  if (thisCut>globallyValid()) { 

3852  // add to global list 

3853  globalCuts_.insert(*thisCut) ; 

3854  } 

3855  } 

3856  for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) { 

3857  whichGenerator_[numberBefore++] = i ; 

3858  const OsiColCut * thisCut = theseCuts.colCutPtr(j) ; 

3859  if (thisCut>globallyValid()) { 

3860  // add to global list 

3861  globalCuts_.insert(*thisCut) ; 

3862  } 

3863  } 

3864  } 

3865  // If at root node and first pass do heuristics without cuts 

3866  if (!numberNodes_&¤tPassNumber_==1) { 

3867  // Save number solutions 

3868  int saveNumberSolutions = numberSolutions_; 

3869  int saveNumberHeuristicSolutions = numberHeuristicSolutions_; 

3870  for (int i = 0;i<numberHeuristics_;i++) { 

3871  // see if heuristic will do anything 

3872  double saveValue = heuristicValue ; 

3873  int ifSol = heuristic_[i]>solution(heuristicValue, 

3874  newSolution); 

3875  if (ifSol>0) { 

3876  // better solution found 

3877  found = i ; 

3878  incrementUsed(newSolution); 

3879  // increment number of solutions so other heuristics can test 

3880  numberSolutions_++; 

3881  numberHeuristicSolutions_++; 

3882  } else { 

3883  heuristicValue = saveValue ; 

3884  } 

3885  } 

3886  // Restore number solutions 

3887  numberSolutions_ = saveNumberSolutions; 

3888  numberHeuristicSolutions_ = saveNumberHeuristicSolutions; 

3889  } 

3890  // Add in any violated saved cuts 

3891  if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) { 

3892  int numberOld = theseCuts.sizeRowCuts(); 

3893  int numberCuts = slackCuts.sizeRowCuts() ; 

3894  int i; 

3895  // possibly extend whichGenerator 

3896  resizeWhichGenerator(numberOld, numberOld+numberCuts); 

3897  for ( i = 0;i<numberCuts;i++) { 

3898  const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ; 

3899  if (thisCut>violated(cbcColSolution_)>100.0*primalTolerance) { 

3900  if (messageHandler()>logLevel()>2) 

3901  printf("Old cut added  violation %g\n", 

3902  thisCut>violated(cbcColSolution_)) ; 

3903  whichGenerator_[numberOld++]=1; 

3904  theseCuts.insert(*thisCut) ; 

3905  } 

3906  } 

3907  } 

3908  /* 

3909  End of the loop to exercise each generator/heuristic. 

3910  

3911  Did any of the heuristics turn up a new solution? Record it before we free 

3912  the vector. 

3913  */ 

3914  if (found >= 0) { 

3915  phase_=4; 

3916  incrementUsed(newSolution); 

3917  setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 

3918  lastHeuristic_ = heuristic_[found]; 

3919  CbcTreeLocal * tree 

3920  = dynamic_cast<CbcTreeLocal *> (tree_); 

3921  if (tree) 

3922  tree>passInSolution(bestSolution_,heuristicValue); 

3923  } 

3924  delete [] newSolution ; 

3925  

3926  #if 0 

3927  // switch on to get all cuts printed 

3928  theseCuts.printCuts() ; 

3929  #endif 

3930  int numberColumnCuts = theseCuts.sizeColCuts() ; 

3931  int numberRowCuts = theseCuts.sizeRowCuts() ; 

3932  if (violated>=0) 

3933  violated = numberRowCuts + numberColumnCuts ; 

3934  /* 

3935  Apply column cuts (aka bound tightening). This may be partially redundant 

3936  for column cuts returned by CglProbing, as generateCuts installs bounds 

3937  from CglProbing when it determines it can fix a variable. 

3938  

3939  TODO: Looks like the use of violated has evolved. The value set above is 

3940  completely ignored. All that's left is violated == 1 indicates some 

3941  cut is violated, violated == 2 indicates infeasibility. Only 

3942  infeasibility warrants exceptional action. 

3943  

3944  TODO: Strikes me that this code will fail to detect infeasibility, because 

3945  the breaks escape the inner loops but the outer loop keeps going. 

3946  Infeasibility in an early cut will be overwritten if a later cut is 

3947  merely violated. 

3948  */ 

3949  if (numberColumnCuts) { 

3950  

3951  #ifdef CBC_DEBUG 

3952  double * oldLower = new double [numberColumns] ; 

3953  double * oldUpper = new double [numberColumns] ; 

3954  memcpy(oldLower,cbcColLower_,numberColumns*sizeof(double)) ; 

3955  memcpy(oldUpper,cbcColUpper_,numberColumns*sizeof(double)) ; 

3956  #endif 

3957  

3958  double integerTolerance = getDblParam(CbcIntegerTolerance) ; 

3959  for (int i = 0;i<numberColumnCuts;i++) { 

3960  const OsiColCut * thisCut = theseCuts.colCutPtr(i) ; 

3961  const CoinPackedVector & lbs = thisCut>lbs() ; 

3962  const CoinPackedVector & ubs = thisCut>ubs() ; 

3963  int j ; 

3964  int n ; 

3965  const int * which ; 

3966  const double * values ; 

3967  n = lbs.getNumElements() ; 

3968  which = lbs.getIndices() ; 

3969  values = lbs.getElements() ; 

3970  for (j = 0;j<n;j++) { 

3971  int iColumn = which[j] ; 

3972  double value = cbcColSolution_[iColumn] ; 

3973  #if CBC_DEBUG>1 

3974  printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn], 

3975  cbcColSolution_[iColumn],oldUpper[iColumn],values[j]) ; 

3976  #endif 

3977  solver_>setColLower(iColumn,values[j]) ; 

3978  if (value<values[j]integerTolerance) 

3979  violated = 1 ; 

3980  if (values[j]>cbcColUpper_[iColumn]+integerTolerance) { 

3981  // infeasible 

3982  violated = 2 ; 

3983  break ; 

3984  } 

3985  } 

3986  n = ubs.getNumElements() ; 

3987  which = ubs.getIndices() ; 

3988  values = ubs.getElements() ; 

3989  for (j = 0;j<n;j++) { 

3990  int iColumn = which[j] ; 

3991  double value = cbcColSolution_[iColumn] ; 

3992  #if CBC_DEBUG>1 

3993  printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn], 

3994  cbcColSolution_[iColumn],oldUpper[iColumn],values[j]) ; 

3995  #endif 

3996  solver_>setColUpper(iColumn,values[j]) ; 

3997  if (value>values[j]+integerTolerance) 

3998  violated = 1 ; 

3999  if (values[j]<cbcColLower_[iColumn]integerTolerance) { 

4000  // infeasible 

4001  violated = 2 ; 

4002  break ; 

4003  } 

4004  } 

4005  } 

4006  #ifdef CBC_DEBUG 

4007  delete [] oldLower ; 

4008  delete [] oldUpper ; 

4009  #endif 

4010  } 

4011  /* 

4012  End installation of column cuts. The break here escapes the numberTries 

4013  loop. 

4014  */ 

4015  if (violated == 2!feasible) { 

4016  // infeasible 

4017  feasible = false ; 

4018  violated = 2; 

4019  if (!numberNodes_) 

4020  messageHandler()>message(CBC_INFEAS, 

4021  messages()) 

4022  << CoinMessageEol ; 

4023  break ; 

4024  } 

4025  /* 

4026  Now apply the row (constraint) cuts. This is a bit more work because we need 

4027  to obtain and augment the current basis. 

4028  

4029  TODO: Why do this work, if there are no row cuts? The current basis will do 

4030  just fine. 

4031  */ 

4032  int numberRowsNow = solver_>getNumRows() ; 

4033  assert(numberRowsNow == numberRowsAtStart+lastNumberCuts) ; 

4034  int numberToAdd = theseCuts.sizeRowCuts() ; 

4035  numberNewCuts_ = lastNumberCuts+numberToAdd ; 

4036  /* 

4037  Now actually add the row cuts and reoptimise. 

4038  

4039  Install the cuts in the solver using applyRowCuts and 

4040  augment the basis with the corresponding slack. We also add each row cut to 

4041  the set of row cuts (cuts.insert()) supplied as a parameter. The new basis 

4042  must be set with setWarmStart(). 

4043  

4044  TODO: Seems to me the original code could allocate addCuts with size 0, if 

4045  numberRowCuts was 0 and numberColumnCuts was nonzero. That might 

4046  explain the memory fault noted in the comment by AJK. Unfortunately, 

4047  just commenting out the delete[] results in massive memory leaks. Try 

4048  a revision to separate the row cut case. Why do we need addCuts at 

4049  all? A typing issue, apparently: OsiCut vs. OsiRowCut. 

4050  

4051  TODO: It looks to me as if numberToAdd and numberRowCuts are identical at 

4052  this point. Confirm & get rid of one of them. 

4053  

4054  TODO: Any reason why the three loops can't be consolidated? 

4055  */ 

4056  if (numberRowCuts > 0  numberColumnCuts > 0) 

4057  { if (numberToAdd > 0) 

4058  { int i ; 

4059  // Faster to add all at once 

4060  const OsiRowCut ** addCuts = new const OsiRowCut * [numberToAdd] ; 

4061  for (i = 0 ; i < numberToAdd ; i++) 

4062  { addCuts[i] = &theseCuts.rowCut(i) ; } 

4063  solver_>applyRowCuts(numberToAdd,addCuts) ; 

4064  // AJK this caused a memory fault on Win32 

4065  // may be okay now with ** form 

4066  delete [] addCuts ; 

4067  for (i = 0 ; i < numberToAdd ; i++) 

4068  { cuts.insert(theseCuts.rowCut(i)) ; } 

4069  CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

4070  assert(basis != NULL); // make sure not volume 

4071  /* dylp bug 

4072  

4073  Consistent size used by OsiDylp as sanity check. Implicit resize seen 

4074  as an error. Hence this call to resize is necessary. 

4075  */ 

4076  basis>resize(numberRowsAtStart+numberNewCuts_,numberColumns) ; 

4077  for (i = 0 ; i < numberToAdd ; i++) 

4078  { basis>setArtifStatus(numberRowsNow+i, 

4079  CoinWarmStartBasis::basic) ; } 

4080  if (solver_>setWarmStart(basis) == false) 

4081  { throw CoinError("Fail setWarmStart() after cut installation.", 

4082  "solveWithCuts","CbcModel") ; } 

4083  delete basis; 

4084  } 

4085  feasible = resolve(node ? node>nodeInfo() : NULL,2) ; 

4086  if ( CoinCpuTime()dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] ) 

4087  numberTries=0; // exit 

4088  # ifdef CBC_DEBUG 

4089  printf("Obj value after cuts %g %d rows\n",solver_>getObjValue(), 

4090  solver_>getNumRows()) ; 

4091  if (onOptimalPath && !solver_>isDualObjectiveLimitReached()) 

4092  assert(feasible) ; 

4093  # endif 

4094  } 

4095  /* 

4096  No cuts. Cut short the cut generation (numberTries) loop. 

4097  */ 

4098  else 

4099  { numberTries = 0 ; } 

4100  /* 

4101  If the problem is still feasible, first, call takeOffCuts() to remove cuts 

4102  that are now slack. takeOffCuts() will call the solver to reoptimise if 

4103  that's needed to restore a valid solution. 

4104  

4105  Next, see if we should quit due to diminishing returns: 

4106  * we've tried three rounds of cut generation and we're getting 

4107  insufficient improvement in the objective; or 

4108  * we generated no cuts; or 

4109  * the solver declared optimality with 0 iterations after we added the 

4110  cuts generated in this round. 

4111  If we decide to keep going, prep for the next iteration. 

4112  

4113  It just seems more safe to tell takeOffCuts() to call resolve(), even if 

4114  we're not continuing cut generation. Otherwise code executed between here 

4115  and final disposition of the node will need to be careful not to access the 

4116  lp solution. It can happen that we lose feasibility in takeOffCuts  

4117  numerical jitters when the cutoff bound is epsilon less than the current 

4118  best, and we're evaluating an alternative optimum. 

4119  

4120  TODO: After successive rounds of code motion, there seems no need to 

4121  distinguish between the two checks for aborting the cut generation 

4122  loop. Confirm and clean up. 

4123  */ 

4124  if (feasible) 

4125  { int cutIterations = solver_>getIterationCount() ; 

4126  if (numberOldActiveCuts_+numberNewCuts_) { 

4127  OsiCuts * saveCuts = node ? NULL : &slackCuts; 

4128  takeOffCuts(cuts,resolveAfterTakeOffCuts_,saveCuts) ; 

4129  if (solver_>isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_) 

4130  { feasible = false ; 

4131  # ifdef CBC_DEBUG 

4132  double z = solver_>getObjValue() ; 

4133  double cut = getCutoff() ; 

4134  printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n", 

4135  zcut,z,cut) ; 

4136  # endif 

4137  } 

4138  } 

4139  if (feasible) { 

4140  numberRowsAtStart = numberOldActiveCuts_+numberRowsAtContinuous_ ; 

4141  lastNumberCuts = numberNewCuts_ ; 

4142  if (direction*solver_>getObjValue() < lastObjective+minimumDrop && 

4143  currentPassNumber_ >= 3) 

4144  { numberTries = 0 ; } 

4145  if (numberRowCuts+numberColumnCuts == 0  cutIterations == 0) 

4146  { break ; } 

4147  if (numberTries > 0) { 

4148  reducedCostFix() ; 

4149  lastObjective = direction*solver_>getObjValue() ; 

4150  } 

4151  } 

4152  } 

4153  /* 

4154  We've lost feasibility  this node won't be referencing the cuts we've 

4155  been collecting, so decrement the reference counts. 

4156  */ 

4157  if (!feasible) 

4158  { int i ; 

4159  for (i = 0;i<currentNumberCuts_;i++) { 

4160  // take off node 

4161  if (addedCuts_[i]) { 

4162  if (!addedCuts_[i]>decrement()) 

4163  delete addedCuts_[i] ; 

4164  addedCuts_[i] = NULL ; 

4165  } 

4166  } 

4167  numberTries = 0 ; 

4168  } 

4169  } while (numberTries>0) ; 

4170  //check feasibility. 

4171  //If solution seems to be integer feasible calling setBestSolution 

4172  //will eventually add extra global cuts which we need to install at 

4173  //the nodes 

4174  

4175  

4176  if(feasible&&solverCharacteristics_>solutionAddsCuts()) //check integer feasibility 

4177  { 

4178  bool integerFeasible = true; 

4179  const double * save = testSolution_; 

4180  testSolution_ = solver_>getColSolution(); 

4181  for (int i=0;i<numberObjects_ && integerFeasible;i++) 

4182  { 

4183  int preferredWay; 

4184  double infeasibility = object_[i]>infeasibility(preferredWay); 

4185  if(infeasibility) 

4186  integerFeasible = false; 

4187  } 

4188  testSolution_ = save; 

4189  if(integerFeasible)//update 

4190  { 

4191  double objValue = solver_>getObjValue(); 

4192  int numberGlobalBefore = globalCuts_.sizeRowCuts(); 

4193  // SOLUTION2 so won't up cutoff or print message 

4194  setBestSolution(CBC_SOLUTION2, objValue, 

4195  solver_>getColSolution(),0); 

4196  int numberGlobalAfter = globalCuts_.sizeRowCuts(); 

4197  int numberToAdd = numberGlobalAfter  numberGlobalBefore; 

4198  if(numberToAdd > 0) 

4199  //We have added some cuts say they are tight at that node 

4200  //Basis and lp should already have been updated 

4201  { 

4202  feasible = (solver_>isProvenOptimal() && 

4203  !solver_>isDualObjectiveLimitReached()) ; 

4204  if(feasible) 

4205  { 

4206  int numberCuts = numberNewCuts_ = cuts.sizeRowCuts(); 

4207  // possibly extend whichGenerator 

4208  resizeWhichGenerator(numberCuts, numberToAdd+numberCuts); 

4209  

4210  for (int i = numberGlobalBefore ; i < numberGlobalAfter ; i++) 

4211  { 

4212  whichGenerator_[numberNewCuts_++]=1; 

4213  cuts.insert(globalCuts_.rowCut(i)) ; 

4214  } 

4215  numberNewCuts_ = lastNumberCuts+numberToAdd; 

4216  //now take off the cuts which are not tight anymore 

4217  takeOffCuts(cuts,resolveAfterTakeOffCuts_, NULL) ; 

4218  if (solver_>isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_) 

4219  { feasible = false ;} 

4220  } 

4221  if(!feasible) //node will be fathomed 

4222  { 

4223  for (int i = 0;i<currentNumberCuts_;i++) 

4224  { 

4225  // take off node 

4226  if (addedCuts_[i]) 

4227  { 

4228  if (!addedCuts_[i]>decrement()) 

4229  delete addedCuts_[i] ; 

4230  addedCuts_[i] = NULL ; 

4231  } 

4232  } 

4233  } 

4234  } 

4235  } 

4236  } 

4237  // Reduced cost fix at end 

4238  if (solver_>isProvenOptimal()) 

4239  reducedCostFix(); 

4240  // If at root node do heuristics 

4241  if (!numberNodes_) { 

4242  // First see if any cuts are slack 

4243  int numberRows = solver_>getNumRows(); 

4244  int numberAdded = numberRowsnumberRowsAtContinuous_; 

4245  if (numberAdded) { 

4246  CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

4247  assert(basis != NULL); 

4248  int * added = new int[numberAdded]; 

4249  int nDelete = 0; 

4250  for (int j=numberRowsAtContinuous_;j<numberRows;j++) { 

4251  if (basis>getArtifStatus(j)==CoinWarmStartBasis::basic) { 

4252  //printf("%d slack!\n",j); 

4253  added[nDelete++]=j; 

4254  } 

4255  } 

4256  if (nDelete) { 

4257  solver_>deleteRows(nDelete,added); 

4258  } 

4259  delete [] added; 

4260  delete basis ; 

4261  } 

4262  // mark so heuristics can tell 

4263  int savePass=currentPassNumber_; 

4264  currentPassNumber_=999999; 

4265  double * newSolution = new double [numberColumns] ; 

4266  double heuristicValue = getCutoff() ; 

4267  int found = 1; // no solution found 

4268  for (int i = 0;i<numberHeuristics_;i++) { 

4269  // see if heuristic will do anything 

4270  double saveValue = heuristicValue ; 

4271  int ifSol = heuristic_[i]>solution(heuristicValue, 

4272  newSolution); 

4273  if (ifSol>0) { 

4274  // better solution found 

4275  found = i ; 

4276  incrementUsed(newSolution); 

4277  } else { 

4278  heuristicValue = saveValue ; 

4279  } 

4280  } 

4281  currentPassNumber_=savePass; 

4282  if (found >= 0) { 

4283  phase_=4; 

4284  incrementUsed(newSolution); 

4285  setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 

4286  lastHeuristic_ = heuristic_[found]; 

4287  } 

4288  delete [] newSolution ; 

4289  } 

4290  // Up change due to cuts 

4291  if (feasible) 

4292  sumChangeObjective2_ += solver_>getObjValue()*solver_>getObjSense() 

4293   objectiveValue ; 

4294  //if ((numberNodes_%100)==0) 

4295  //printf("XXb sum obj changed by %g\n",sumChangeObjective2_); 

4296  /* 

4297  End of cut generation loop. 

4298  

4299  Now, consider if we want to disable or adjust the frequency of use for any 

4300  of the cut generators. If the client specified a positive number for 

4301  howOften, it will never change. If the original value was negative, it'll 

4302  be converted to 1000000+howOften, and this value will be adjusted each 

4303  time fullScan is true. Actual cut generation is performed every 

4304  howOften%1000000 nodes; the 1000000 offset is just a convenient way to 

4305  specify that the frequency is adjustable. 

4306  

4307  During cut generation, we recorded the number of cuts produced by each 

4308  generator for this node. For all cuts, whichGenerator records the generator 

4309  that produced a cut. 

4310  

4311  TODO: All this should probably be hidden in a method of the CbcCutGenerator 

4312  class. 

4313  

4314  TODO: Can the loop that scans over whichGenerator to accumulate per generator 

4315  counts be replaced by values in countRowCuts and countColumnCuts? 

4316  */ 

4317  #ifdef NODE_LOG 

4318  int fatherNum = (node == NULL) ? 1 : node>nodeNumber(); 

4319  double value = (node == NULL) ? 1 : node>branchingObject()>value(); 

4320  string bigOne = (solver_>getIterationCount() > 30) ? "*******" : ""; 

4321  string way = (node == NULL) ? "" : (node>branchingObject()>way())==1 ? "Down" : "Up"; 

4322  std::cout<<"Node "<<numberNodes_<<", father "<<fatherNum<<", #iterations "<<solver_>getIterationCount()<<", sol value : "<<solver_>getObjValue()<<std::endl; 

4323  #endif 

4324  if (fullScan&&numberCutGenerators_) { 

4325  /* If cuts just at root node then it will probably be faster to 

4326  update matrix and leave all in */ 

4327  int willBeCutsInTree=0; 

4328  double thisObjective = solver_>getObjValue()*direction ; 

4329  // get sizes 

4330  int numberRowsAdded = solver_>getNumRows()numberRowsAtStart; 

4331  CoinBigIndex numberElementsAdded = solver_>getNumElements()numberElementsAtStart ; 

4332  double densityOld = ((double) numberElementsAtStart)/((double) numberRowsAtStart); 

4333  double densityNew = numberRowsAdded ? ((double) (numberElementsAdded))/((double) numberRowsAdded) 

4334  : 0.0; 

4335  if (!numberNodes_) { 

4336  if (numberRowsAdded) 

4337  handler_>message(CBC_CUTS_STATS,messages_) 

4338  <<numberRowsAdded 

4339  <<densityNew 

4340  <<CoinMessageEol ; 

4341  if (thisObjectivestartObjective<1.0e5&&numberElementsAdded>0.2*numberElementsAtStart) 

4342  willBeCutsInTree=1; 

4343  } 

4344  if ((numberRowsAdded>100+0.5*numberRowsAtStart 

4345  numberElementsAdded>0.5*numberElementsAtStart) 

4346  &&(densityNew>200.0&&numberRowsAdded>100&&densityNew>2.0*densityOld)) { 

4347  // much bigger 

4348  //if (thisObjectivestartObjective<0.1*fabs(startObjective)+1.0e5) 

4349  //willBeCutsInTree=1; 

4350  //printf("Cuts will be taken off , %d rows added with density %g\n", 

4351  // numberRowsAdded,densityNew); 

4352  } 

4353  if (densityNew>100.0&&numberRowsAdded>2&&densityNew>2.0*densityOld) { 

4354  //if (thisObjectivestartObjective<0.1*fabs(startObjective)+1.0e5) 

4355  //willBeCutsInTree=2; 

4356  //printf("Density says no cuts ? , %d rows added with density %g\n", 

4357  // numberRowsAdded,densityNew); 

4358  } 

4359  // Root node or every so often  see what to turn off 

4360  int i ; 

4361  for (i = 0;i<numberCutGenerators_;i++) { 

4362  int howOften = generator_[i]>howOften() ; 

4363  if (howOften>90) 

4364  willBeCutsInTree=0; 

4365  } 

4366  if (!numberNodes_) 

4367  handler_>message(CBC_ROOT,messages_) 

4368  <<numberNewCuts_ 

4369  <<startObjective<<thisObjective 

4370  <<currentPassNumber_ 

4371  <<CoinMessageEol ; 

4372  int * count = new int[numberCutGenerators_] ; 

4373  memset(count,0,numberCutGenerators_*sizeof(int)) ; 

4374  int numberActiveGenerators=0; 

4375  for (i = 0;i<numberNewCuts_;i++) { 

4376  int iGenerator = whichGenerator_[i]; 

4377  if (iGenerator>=0&&iGenerator<numberCutGenerators_) 

4378  count[iGenerator]++ ; 

4379  } 

4380  double totalCuts = 0.0 ; 

4381  //#define JUST_ACTIVE 

4382  for (i = 0;i<numberCutGenerators_;i++) { 

4383  if (countRowCuts[i]countColumnCuts[i]) 

4384  numberActiveGenerators++; 

4385  #ifdef JUST_ACTIVE 

4386  totalCuts += count[i] + 5.0*countColumnCuts[i] ; 

4387  #else 

4388  totalCuts += countRowCuts[i] + 5.0*countColumnCuts[i] ; 

4389  #endif 

4390  } 

4391  double small = (0.5* totalCuts) / 

4392  ((double) numberActiveGenerators) ; 

4393  for (i = 0;i<numberCutGenerators_;i++) { 

4394  int howOften = generator_[i]>howOften() ; 

4395  if (willBeCutsInTree<0&&howOften==98) 

4396  howOften =99; 

4397  if (howOften==98&&generator_[i]>switchOffIfLessThan()<0) { 

4398  if (thisObjectivestartObjective<0.005*fabs(startObjective)+1.0e5) 

4399  howOften=99; // switch off 

4400  if (thisObjectivestartObjective<0.1*fabs(startObjective)+1.0e5 

4401  &&5*solver_>getNumRows()<solver_>getNumCols()) 

4402  howOften=99; // switch off 

4403  } 

4404  if (howOften<99) 

4405  continue ; 

4406  if (howOften<0howOften >= 1000000) { 

4407  if( !numberNodes_) { 

4408  // If small number switch mostly off 

4409  #ifdef JUST_ACTIVE 

4410  double thisCuts = count[i] + 5.0*countColumnCuts[i] ; 

4411  #else 

4412  double thisCuts = countRowCuts[i] + 5.0*countColumnCuts[i] ; 

4413  #endif 

4414  if (!thisCutshowOften == 99) { 

4415  if (howOften == 99howOften == 98) 

4416  howOften = 100 ; 

4417  else 

4418  howOften = 1000000+SCANCUTS; // wait until next time 

4419  } else if (thisCuts<small) { 

4420  int k = (int) sqrt(small/thisCuts) ; 

4421  if (howOften!=98) 

4422  howOften = k+1000000 ; 

4423  else 

4424  howOften=100; 

4425  } else { 

4426  howOften = 1+1000000 ; 

4427  if (thisObjectivestartObjective<0.1*fabs(startObjective)+1.0e5) 

4428  generator_[i]>setWhatDepth(5); 

4429  } 

4430  } 

4431  // If cuts useless switch off 

4432  if (numberNodes_>=100000&&sumChangeObjective1_>2.0e2*(sumChangeObjective2_+1.0e12)) { 

4433  howOften = 1000000+SCANCUTS; // wait until next time 

4434  //printf("switch off cut %d due to lack of use\n",i); 

4435  } 

4436  } 

4437  if (howOften>=0&&generator_[i]>generator()>mayGenerateRowCutsInTree()) 

4438  willBeCutsInTree=1; 

4439  

4440  generator_[i]>setHowOften(howOften) ; 

4441  if (howOften>=1000000&&howOften<2000000&&0) { 

4442  // Go to depth 

4443  int bias=1; 

4444  if (howOften==1+1000000) 

4445  generator_[i]>setWhatDepth(bias+1); 

4446  else if (howOften<=10+1000000) 

4447  generator_[i]>setWhatDepth(bias+2); 

4448  else 

4449  generator_[i]>setWhatDepth(bias+1000); 

4450  } 

4451  int newFrequency = generator_[i]>howOften()%1000000 ; 

4452  // increment cut counts 

4453  generator_[i]>incrementNumberCutsInTotal(countRowCuts[i]); 

4454  generator_[i]>incrementNumberCutsActive(count[i]); 

4455  if (handler_>logLevel()>1!numberNodes_) { 

4456  handler_>message(CBC_GENERATOR,messages_) 

4457  <<i 

4458  <<generator_[i]>cutGeneratorName() 

4459  <<countRowCuts[i]<<count[i] 

4460  <<countColumnCuts[i]; 

4461  handler_>printing(!numberNodes_&&generator_[i]>timing()) 

4462  <<generator_[i]>timeInCutGenerator(); 

4463  handler_>message() 

4464  <<newFrequency 

4465  <<CoinMessageEol ; 

4466  } 

4467  } 

4468  delete [] count ; 

4469  if( !numberNodes_) { 

4470  if (willBeCutsInTree==2) 

4471  willBeCutsInTree=0; 

4472  if( willBeCutsInTree<=0) { 

4473  // Take off cuts 

4474  cuts = OsiCuts(); 

4475  numberNewCuts_=0; 

4476  if (!willBeCutsInTree) { 

4477  // update size of problem 

4478  numberRowsAtContinuous_ = solver_>getNumRows() ; 

4479  } else { 

4480  // take off cuts 

4481  int numberRows = solver_>getNumRows(); 

4482  int numberAdded = numberRowsnumberRowsAtContinuous_; 

4483  if (numberAdded) { 

4484  int * added = new int[numberAdded]; 

4485  for (int i=0;i<numberAdded;i++) 

4486  added[i]=i+numberRowsAtContinuous_; 

4487  solver_>deleteRows(numberAdded,added); 

4488  delete [] added; 

4489  // resolve so optimal 

4490  solver_>resolve(); 

4491  } 

4492  } 

4493  #ifdef CBC_USE_CLP 

4494  OsiClpSolverInterface * clpSolver 

4495  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

4496  if (clpSolver) { 

4497  // Maybe solver might like to know only column bounds will change 

4498  //int options = clpSolver>specialOptions(); 

4499  //clpSolver>setSpecialOptions(options128); 

4500  clpSolver>synchronizeModel(); 

4501  } 

4502  #endif 

4503  } else { 

4504  #ifdef CBC_USE_CLP 

4505  OsiClpSolverInterface * clpSolver 

4506  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

4507  if (clpSolver) { 

4508  // make sure factorization can't carry over 

4509  int options = clpSolver>specialOptions(); 

4510  clpSolver>setSpecialOptions(options&(~8)); 

4511  } 

4512  #endif 

4513  } 

4514  } 

4515  } else if (numberCutGenerators_) { 

4516  int i; 

4517  // add to counts anyway 

4518  for (i = 0;i<numberCutGenerators_;i++) 

4519  generator_[i]>incrementNumberCutsInTotal(countRowCuts[i]); 

4520  // What if not feasible as cuts may have helped 

4521  if (feasible) { 

4522  for (i = 0;i<numberNewCuts_;i++) { 

4523  int iGenerator = whichGenerator_[i]; 

4524  if (iGenerator>=0) 

4525  generator_[iGenerator]>incrementNumberCutsActive(); 

4526  } 

4527  } 

4528  } 

4529  

4530  delete [] countRowCuts ; 

4531  delete [] countColumnCuts ; 

4532  

4533  

4534  #ifdef CHECK_CUT_COUNTS 

4535  if (feasible) 

4536  { 

4537  CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

4538  printf("solveWithCuts: Number of rows at end (only active cuts) %d\n", 

4539  numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ; 

4540  basis>print() ; delete basis;} 

4541  #endif 

4542  #ifdef CBC_DEBUG 

4543  if (onOptimalPath && !solver_>isDualObjectiveLimitReached()) 

4544  assert(feasible) ; 

4545  #endif 

4546  

4547  return feasible ; } 

4548  

4549  

4550  /* 

4551  Remove slack cuts. We obtain a basis and scan it. Cuts with basic slacks 

4552  are purged. If any cuts are purged, resolve() is called to restore the 

4553  solution held in the solver. If resolve() pivots, there's the possibility 

4554  that a slack may be pivoted in (trust me :), so the process iterates. 

4555  Setting allowResolve to false will suppress reoptimisation (but see note 

4556  below). 

4557  

4558  At the level of the solver's constraint system, loose cuts are really 

4559  deleted. There's an implicit assumption that deleteRows will also update 

4560  the active basis in the solver. 

4561  

4562  At the level of nodes and models, it's more complicated. 

4563  

4564  New cuts exist only in the collection of cuts passed as a parameter. They 

4565  are deleted from the collection and that's the end of them. 

4566  

4567  Older cuts have made it into addedCuts_. Two separate actions are needed. 

4568  The reference count for the CbcCountRowCut object is decremented. If this 

4569  count falls to 0, the node which owns the cut is located, the reference to 

4570  the cut is removed, and then the cut object is destroyed (courtesy of the 

4571  CbcCountRowCut destructor). We also need to set the addedCuts_ entry to 

4572  NULL. This is important so that when it comes time to generate basis edits 

4573  we can tell this cut was dropped from the basis during processing of the 

4574  node. 

4575  

4576  NOTE: In general, it's necessary to call resolve() after purging slack 

4577  cuts. Deleting constraints constitutes a change in the problem, and 

4578  an OSI is not required to maintain a valid solution when the problem 

4579  is changed. But ... it's really useful to maintain the active basis, 

4580  and the OSI is supposed to do that. (Yes, it's splitting hairs.) In 

4581  some places, it's possible to know that the solution will never be 

4582  consulted after this call, only the basis. (E.g., this routine is 

4583  called as a last act before generating info to place the node in the 

4584  live set.) For such use, set allowResolve to false. 

4585  

4586  TODO: No real harm would be done if we just ignored the rare occasion when 

4587  the call to resolve() pivoted a slack back into the basis. It's a 

4588  minor inefficiency, at worst. But it does break assertions which 

4589  check that there are no loose cuts in the basis. It might be better 

4590  to remove the assertions. 

4591  */ 

4592  

4593  void 

4594  CbcModel::takeOffCuts (OsiCuts &newCuts, 

4595  bool allowResolve, OsiCuts * saveCuts) 

4596  

4597  { // int resolveIterations = 0 ; 

4598  int firstOldCut = numberRowsAtContinuous_ ; 

4599  int totalNumberCuts = numberNewCuts_+numberOldActiveCuts_ ; 

4600  int *solverCutIndices = new int[totalNumberCuts] ; 

4601  int *newCutIndices = new int[numberNewCuts_] ; 

4602  const CoinWarmStartBasis* ws ; 

4603  CoinWarmStartBasis::Status status ; 

4604  bool needPurge = true ; 

4605  /* 

4606  The outer loop allows repetition of purge in the event that reoptimisation 

4607  changes the basis. To start an iteration, clear the deletion counts and grab 

4608  the current basis. 

4609  */ 

4610  while (needPurge) 

4611  { int numberNewToDelete = 0 ; 

4612  int numberOldToDelete = 0 ; 

4613  int i ; 

4614  ws = dynamic_cast<const CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

4615  /* 

4616  Scan the basis entries of the old cuts generated prior to this round of cut 

4617  generation. Loose cuts are `removed' by decrementing their reference count 

4618  and setting the addedCuts_ entry to NULL. (If the reference count falls to 

4619  0, they're really deleted. See CbcModel and CbcCountRowCut doc'n for 

4620  principles of cut handling.) 

4621  */ 

4622  int oldCutIndex = 0 ; 

4623  for (i = 0 ; i < numberOldActiveCuts_ ; i++) 

4624  { status = ws>getArtifStatus(i+firstOldCut) ; 

4625  while (!addedCuts_[oldCutIndex]) oldCutIndex++ ; 

4626  assert(oldCutIndex < currentNumberCuts_) ; 

4627  // always leave if from nextRowCut_ 

4628  if (status == CoinWarmStartBasis::basic&& 

4629  addedCuts_[oldCutIndex]>effectiveness()!=COIN_DBL_MAX) 

4630  { solverCutIndices[numberOldToDelete++] = i+firstOldCut ; 

4631  if (saveCuts) { 

4632  // send to cut pool 

4633  OsiRowCut * slackCut = addedCuts_[oldCutIndex]; 

4634  if (slackCut>effectiveness()!=1.234) { 

4635  slackCut>setEffectiveness(1.234); 

4636  saveCuts>insert(*slackCut); 

4637  } 

4638  } 

4639  if (addedCuts_[oldCutIndex]>decrement() == 0) 

4640  delete addedCuts_[oldCutIndex] ; 

4641  addedCuts_[oldCutIndex] = NULL ; 

4642  oldCutIndex++ ; } 

4643  else 

4644  { oldCutIndex++ ; } } 

4645  /* 

4646  Scan the basis entries of the new cuts generated with this round of cut 

4647  generation. At this point, newCuts is the only record of the new cuts, so 

4648  when we delete loose cuts from newCuts, they're really gone. newCuts is a 

4649  vector, so it's most efficient to compress it (eraseRowCut) from back to 

4650  front. 

4651  */ 

4652  int firstNewCut = firstOldCut+numberOldActiveCuts_ ; 

4653  int k = 0 ; 

4654  for (i = 0 ; i < numberNewCuts_ ; i++) 

4655  { status = ws>getArtifStatus(i+firstNewCut) ; 

4656  if (status == CoinWarmStartBasis::basic&&whichGenerator_[i]!=2) 

4657  { solverCutIndices[numberNewToDelete+numberOldToDelete] = i+firstNewCut ; 

4658  newCutIndices[numberNewToDelete++] = i ; } 

4659  else 

4660  { // save which generator did it 

4661  whichGenerator_[k++] = whichGenerator_[i] ; } } 

4662  delete ws ; 

4663  for (i = numberNewToDelete1 ; i >= 0 ; i) 

4664  { int iCut = newCutIndices[i] ; 

4665  if (saveCuts) { 

4666  // send to cut pool 

4667  OsiRowCut * slackCut = newCuts.rowCutPtr(iCut); 

4668  if (slackCut>effectiveness()!=1.234) { 

4669  slackCut>setEffectiveness(1.234); 

4670  saveCuts>insert(*slackCut); 

4671  } 

4672  } 

4673  newCuts.eraseRowCut(iCut) ; } 

4674  /* 

4675  Did we delete anything? If so, delete the cuts from the constraint system 

4676  held in the solver and reoptimise unless we're forbidden to do so. If the 

4677  call to resolve() results in pivots, there's the possibility we again have 

4678  basic slacks. Repeat the purging loop. 

4679  */ 

4680  if (numberNewToDelete+numberOldToDelete > 0) 

4681  { solver_>deleteRows(numberNewToDelete+numberOldToDelete, 

4682  solverCutIndices) ; 

4683  numberNewCuts_ = numberNewToDelete ; 

4684  numberOldActiveCuts_ = numberOldToDelete ; 

4685  # ifdef CBC_DEBUG 

4686  printf("takeOffCuts: purged %d+%d cuts\n", numberOldToDelete, 

4687  numberNewToDelete ); 

4688  # endif 

4689  if (allowResolve) 

4690  { 

4691  phase_=3; 

4692  // can do quick optimality check 

4693  int easy=2; 

4694  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ; 

4695  solver_>resolve() ; 

4696  setPointers(solver_); 

4697  solver_>setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ; 

4698  if (solver_>getIterationCount() == 0) 

4699  { needPurge = false ; } 

4700  # ifdef CBC_DEBUG 

4701  else 

4702  { printf( "Repeating purging loop. %d iters.\n", 

4703  solver_>getIterationCount()); } 

4704  # endif 

4705  } 

4706  else 

4707  { needPurge = false ; } } 

4708  else 

4709  { needPurge = false ; } } 

4710  /* 

4711  Clean up and return. 

4712  */ 

4713  delete [] solverCutIndices ; 

4714  delete [] newCutIndices ; 

4715  } 

4716  

4717  

4718  

4719  int 

4720  CbcModel::resolve(CbcNodeInfo * parent, int whereFrom) 

4721  { 

4722  // We may have deliberately added in violated cuts  check to avoid message 

4723  int iRow; 

4724  int numberRows = solver_>getNumRows(); 

4725  const double * rowLower = solver_>getRowLower(); 

4726  const double * rowUpper = solver_>getRowUpper(); 

4727  bool feasible=true; 

4728  for (iRow= numberRowsAtContinuous_;iRow<numberRows;iRow++) { 

4729  if (rowLower[iRow]>rowUpper[iRow]+1.0e8) 

4730  feasible=false; 

4731  } 

4732  // Can't happen if strong branching as would have been found before 

4733  if (!numberStrong_&&numberObjects_>numberIntegers_) { 

4734  int iColumn; 

4735  int numberColumns = solver_>getNumCols(); 

4736  const double * columnLower = solver_>getColLower(); 

4737  const double * columnUpper = solver_>getColUpper(); 

4738  for (iColumn= 0;iColumn<numberColumns;iColumn++) { 

4739  if (columnLower[iColumn]>columnUpper[iColumn]+1.0e5) 

4740  feasible=false; 

4741  } 

4742  } 

4743  #if 0 

4744  { 

4745  int iColumn; 

4746  int numberColumns = solver_>getNumCols(); 

4747  const double * columnLower = solver_>getColLower(); 

4748  const double * columnUpper = solver_>getColUpper(); 

4749  const double * sol = solver_>getColSolution(); 

4750  bool feasible=true; 

4751  double xx[37]; 

4752  memset(xx,0,sizeof(xx)); 

4753  xx[6]=xx[7]=xx[9]=xx[18]=xx[26]=xx[36]=1.0; 

4754  for (iColumn= 0;iColumn<numberColumns;iColumn++) { 

4755  if (solver_>isInteger(iColumn)) { 

4756  //printf("yy %d at %g (bounds %g, %g)\n",iColumn, 

4757  // sol[iColumn],columnLower[iColumn],columnUpper[iColumn]); 

4758  if (columnLower[iColumn]>xx[iColumn] 

4759  columnUpper[iColumn]<xx[iColumn]) 

4760  feasible=false; 

4761  //solver_>setColLower(iColumn,CoinMin(xx[iColumn],columnUpper[iColumn])); 

4762  //solver_>setColUpper(iColumn,CoinMax(xx[iColumn],columnLower[iColumn])); 

4763  } 

4764  } 

4765  if (feasible) 

4766  printf("On Feasible path\n"); 

4767  } 

4768  #endif 

4769  /* 

4770  Reoptimize. Consider the possibility that we should fathom on bounds. But be 

4771  careful  where the objective takes on integral values, we may want to keep 

4772  a solution where the objective is right on the cutoff. 

4773  */ 

4774  if (feasible) 

4775  { 

4776  solver_>resolve() ; 

4777  numberIterations_ += solver_>getIterationCount() ; 

4778  feasible = (solver_>isProvenOptimal() && 

4779  !solver_>isDualObjectiveLimitReached()) ; 

4780  } 

4781  if (0&&feasible) { 

4782  const double * lb = solver_>getColLower(); 

4783  const double * ub = solver_>getColUpper(); 

4784  const double * x = solver_>getColSolution(); 

4785  const double * dj = solver_>getReducedCost(); 

4786  int numberColumns = solver_>getNumCols(); 

4787  for (int i=0;i<numberColumns;i++) { 

4788  if (dj[i]>1.0e4&&ub[i]lb[i]>1.0e4&&x[i]>lb[i]+1.0e4) 

4789  printf("error %d %g %g %g %g\n",i,dj[i],lb[i],x[i],ub[i]); 

4790  if (dj[i]<1.0e4&&ub[i]lb[i]>1.0e4&&x[i]<ub[i]1.0e4) 

4791  printf("error %d %g %g %g %g\n",i,dj[i],lb[i],x[i],ub[i]); 

4792  } 

4793  } 

4794  if (!feasible&& continuousObjective_ <1.0e30) { 

4795  // at root node  double double check 

4796  bool saveTakeHint; 

4797  OsiHintStrength saveStrength; 

4798  solver_>getHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength); 

4799  if (saveTakeHintsaveStrength==OsiHintIgnore) { 

4800  solver_>setHintParam(OsiDoDualInResolve,false,OsiHintDo) ; 

4801  solver_>resolve(); 

4802  solver_>setHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength); 

4803  numberIterations_ += solver_>getIterationCount() ; 

4804  feasible = solver_>isProvenOptimal(); 

4805  // solver_>writeMps("infeas"); 

4806  } 

4807  } 

4808  setPointers(solver_); 

4809  int returnStatus = feasible ? 1 : 0; 

4810  if (strategy_) { 

4811  // user can play clever tricks here 

4812  int status = strategy_>status(this,parent,whereFrom); 

4813  if (status>=0) { 

4814  if (status==0) 

4815  returnStatus = 1; 

4816  else if(status==1) 

4817  returnStatus=1; 

4818  else 

4819  returnStatus=0; 

4820  } 

4821  } 

4822  return returnStatus ; 

4823  } 

4824  

4825  

4826  /* Set up objects. Only do ones whose length is in range. 

4827  If makeEquality true then a new model may be returned if 

4828  modifications had to be made, otherwise "this" is returned. 

4829  

4830  Could use Probing at continuous to extend objects 

4831  */ 

4832  CbcModel * 

4833  CbcModel::findCliques(bool makeEquality, 

4834  int atLeastThisMany, int lessThanThis, int defaultValue) 

4835  { 

4836  // No objects are allowed to exist 

4837  assert(numberObjects_==numberIntegers_!numberObjects_); 

4838  CoinPackedMatrix matrixByRow(*solver_>getMatrixByRow()); 

4839  int numberRows = solver_>getNumRows(); 

4840  int numberColumns = solver_>getNumCols(); 

4841  

4842  // We may want to add columns 

4843  int numberSlacks=0; 

4844  int * rows = new int[numberRows]; 

4845  double * element =new double[numberRows]; 

4846  

4847  int iRow; 

4848  

4849  findIntegers(true); 

4850  numberObjects_=numberIntegers_; 

4851  

4852  int numberCliques=0; 

4853  CbcObject ** object = new CbcObject * [numberRows]; 

4854  int * which = new int[numberIntegers_]; 

4855  char * type = new char[numberIntegers_]; 

4856  int * lookup = new int[numberColumns]; 

4857  int i; 

4858  for (i=0;i<numberColumns;i++) 

4859  lookup[i]=1; 

4860  for (i=0;i<numberIntegers_;i++) 

4861  lookup[integerVariable_[i]]=i; 

4862  

4863  // Statistics 

4864  int totalP1=0,totalM1=0; 

4865  int numberBig=0,totalBig=0; 

4866  int numberFixed=0; 

4867  

4868  // Row copy 

4869  const double * elementByRow = matrixByRow.getElements(); 

4870  const int * column = matrixByRow.getIndices(); 

4871  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts(); 

4872  const int * rowLength = matrixByRow.getVectorLengths(); 

4873  

4874  // Column lengths for slacks 

4875  const int * columnLength = solver_>getMatrixByCol()>getVectorLengths(); 

4876  

4877  const double * lower = getColLower(); 

4878  const double * upper = getColUpper(); 

4879  const double * rowLower = getRowLower(); 

4880  const double * rowUpper = getRowUpper(); 

4881  

4882  for (iRow=0;iRow<numberRows;iRow++) { 

4883  int numberP1=0, numberM1=0; 

4884  int j; 

4885  double upperValue=rowUpper[iRow]; 

4886  double lowerValue=rowLower[iRow]; 

4887  bool good=true; 

4888  int slack = 1; 

4889  for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) { 

4890  int iColumn = column[j]; 

4891  int iInteger=lookup[iColumn]; 

4892  if (upper[iColumn]lower[iColumn]<1.0e8) { 

4893  // fixed 

4894  upperValue = lower[iColumn]*elementByRow[j]; 

4895  lowerValue = lower[iColumn]*elementByRow[j]; 

4896  continue; 

4897  } else if (upper[iColumn]!=1.0lower[iColumn]!=0.0) { 

4898  good = false; 

4899  break; 

4900  } else if (iInteger<0) { 

4901  good = false; 

4902  break; 

4903  } else { 

4904  if (columnLength[iColumn]==1) 

4905  slack = iInteger; 

4906  } 

4907  if (fabs(elementByRow[j])!=1.0) { 

4908  good=false; 

4909  break; 

4910  } else if (elementByRow[j]>0.0) { 

4911  which[numberP1++]=iInteger; 

4912  } else { 

4913  numberM1++; 

4914  which[numberIntegers_numberM1]=iInteger; 

4915  } 

4916  } 

4917  int iUpper = (int) floor(upperValue+1.0e5); 

4918  int iLower = (int) ceil(lowerValue1.0e5); 

4919  int state=0; 

4920  if (upperValue<1.0e6) { 

4921  if (iUpper==1numberM1) 

4922  state=1; 

4923  else if (iUpper==numberM1) 

4924  state=2; 

4925  else if (iUpper<numberM1) 

4926  state=3; 

4927  } 

4928  if (!state&&lowerValue>1.0e6) { 

4929  if (iLower==1numberP1) 

4930  state=1; 

4931  else if (iLower==numberP1) 

4932  state=2; 

4933  else if (iLower<numberP1) 

4934  state=3; 

4935  } 

4936  if (good&&state) { 

4937  if (abs(state)==3) { 

4938  // infeasible 

4939  numberObjects_=1; 

4940  break; 

4941  } else if (abs(state)==2) { 

4942  // we can fix all 

4943  numberFixed += numberP1+numberM1; 

4944  if (state>0) { 

4945  // fix all +1 at 0, 1 at 1 

4946  for (i=0;i<numberP1;i++) 

4947  solver_>setColUpper(integerVariable_[which[i]],0.0); 

4948  for (i=0;i<numberM1;i++) 

4949  solver_>setColLower(integerVariable_[which[numberIntegers_i1]], 

4950  1.0); 

4951  } else { 

4952  // fix all +1 at 1, 1 at 0 

4953  for (i=0;i<numberP1;i++) 

4954  solver_>setColLower(integerVariable_[which[i]],1.0); 

4955  for (i=0;i<numberM1;i++) 

4956  solver_>setColUpper(integerVariable_[which[numberIntegers_i1]], 

4957  0.0); 

4958  } 

4959  } else { 

4960  int length = numberP1+numberM1; 

4961  if (length >= atLeastThisMany&&length<lessThanThis) { 

4962  // create object 

4963  bool addOne=false; 

4964  int objectType; 

4965  if (iLower==iUpper) { 

4966  objectType=1; 

4967  } else { 

4968  if (makeEquality) { 

4969  objectType=1; 

4970  element[numberSlacks]=state; 

4971  rows[numberSlacks++]=iRow; 

4972  addOne=true; 

4973  } else { 

4974  objectType=0; 

4975  } 

4976  } 

4977  if (state>0) { 

4978  totalP1 += numberP1; 

4979  totalM1 += numberM1; 

4980  for (i=0;i<numberP1;i++) 

4981  type[i]=1; 

4982  for (i=0;i<numberM1;i++) { 

4983  which[numberP1]=which[numberIntegers_i1]; 

4984  type[numberP1++]=0; 

4985  } 

4986  } else { 

4987  totalP1 += numberM1; 

4988  totalM1 += numberP1; 

4989  for (i=0;i<numberP1;i++) 

4990  type[i]=0; 

4991  for (i=0;i<numberM1;i++) { 

4992  which[numberP1]=which[numberIntegers_i1]; 

4993  type[numberP1++]=1; 

4994  } 

4995  } 

4996  if (addOne) { 

4997  // add in slack 

4998  which[numberP1]=numberIntegers_+numberSlacks1; 

4999  slack = numberP1; 

5000  type[numberP1++]=1; 

5001  } else if (slack >= 0) { 

5002  for (i=0;i<numberP1;i++) { 

5003  if (which[i]==slack) { 

5004  slack=i; 

5005  } 

5006  } 

5007  } 

5008  object[numberCliques] = new CbcClique(this,objectType,numberP1, 

5009  which,type, 

5010  1000000+numberCliques,slack); 

5011  numberCliques++; 

5012  } else if (numberP1+numberM1 >= lessThanThis) { 

5013  // too big 

5014  numberBig++; 

5015  totalBig += numberP1+numberM1; 

5016  } 

5017  } 

5018  } 

5019  } 

5020  delete [] which; 

5021  delete [] type; 

5022  delete [] lookup; 

5023  if (numberCliques<0) { 

5024  printf("*** Problem infeasible\n"); 

5025  } else { 

5026  if (numberCliques) 

5027  printf("%d cliques of average size %g found, %d P1, %d M1\n", 

5028  numberCliques, 

5029  ((double)(totalP1+totalM1))/((double) numberCliques), 

5030  totalP1,totalM1); 

5031  else 

5032  printf("No cliques found\n"); 

5033  if (numberBig) 

5034  printf("%d large cliques ( >= %d) found, total %d\n", 

5035  numberBig,lessThanThis,totalBig); 

5036  if (numberFixed) 

5037  printf("%d variables fixed\n",numberFixed); 

5038  } 

5039  if (numberCliques>0&&numberSlacks&&makeEquality) { 

5040  printf("adding %d integer slacks\n",numberSlacks); 

5041  // add variables to make equality rows 

5042  int * temp = new int[numberIntegers_+numberSlacks]; 

5043  memcpy(temp,integerVariable_,numberIntegers_*sizeof(int)); 

5044  // Get new model 

5045  CbcModel * newModel = new CbcModel(*this); 

5046  OsiSolverInterface * newSolver = newModel>solver(); 

5047  for (i=0;i<numberSlacks;i++) { 

5048  temp[i+numberIntegers_]=i+numberColumns; 

5049  int iRow = rows[i]; 

5050  double value = element[i]; 

5051  double lowerValue = 0.0; 

5052  double upperValue = 1.0; 

5053  double objValue = 0.0; 

5054  CoinPackedVector column(1,&iRow,&value); 

5055  newSolver>addCol(column,lowerValue,upperValue,objValue); 

5056  // set integer 

5057  newSolver>setInteger(numberColumns+i); 

5058  if (value >0) 

5059  newSolver>setRowLower(iRow,rowUpper[iRow]); 

5060  else 

5061  newSolver>setRowUpper(iRow,rowLower[iRow]); 

5062  } 

5063  // replace list of integers 

5064  for (i=0;i<newModel>numberObjects_;i++) 

5065  delete newModel>object_[i]; 

5066  newModel>numberObjects_ = 0; 

5067  delete [] newModel>object_; 

5068  newModel>object_=NULL; 

5069  newModel>findIntegers(true); //Set up all integer objects 

5070  for (i=0;i<numberIntegers_;i++) { 

5071  newModel>modifiableObject(i)>setPriority(object_[i]>priority()); 

5072  } 

5073  if (originalColumns_) { 

5074  // old model had originalColumns 

5075  delete [] newModel>originalColumns_; 

5076  newModel>originalColumns_ = new int[numberColumns+numberSlacks]; 

5077  memcpy(newModel>originalColumns_,originalColumns_,numberColumns*sizeof(int)); 

5078  // mark as not in previous model 

5079  for (i=numberColumns;i<numberColumns+numberSlacks;i++) 

5080  newModel>originalColumns_[i]=1; 

5081  } 

5082  delete [] rows; 

5083  delete [] element; 

5084  newModel>addObjects(numberCliques,object); 

5085  for (;i<numberCliques;i++) 

5086  delete object[i]; 

5087  delete [] object; 

5088  newModel>synchronizeModel(); 

5089  return newModel; 

5090  } else { 

5091  if (numberCliques>0) { 

5092  addObjects(numberCliques,object); 

5093  for (;i<numberCliques;i++) 

5094  delete object[i]; 

5095  synchronizeModel(); 

5096  } 

5097  delete [] object; 

5098  delete [] rows; 

5099  delete [] element; 

5100  return this; 

5101  } 

5102  } 

5103  // Fill in useful estimates 

5104  void 

5105  CbcModel::pseudoShadow(double * down, double * up) 

5106  { 

5107  // Column copy of matrix 

5108  const double * element = solver_>getMatrixByCol()>getElements(); 

5109  const int * row = solver_>getMatrixByCol()>getIndices(); 

5110  const CoinBigIndex * columnStart = solver_>getMatrixByCol()>getVectorStarts(); 

5111  const int * columnLength = solver_>getMatrixByCol()>getVectorLengths(); 

5112  const double *objective = solver_>getObjCoefficients() ; 

5113  int numberColumns = solver_>getNumCols() ; 

5114  double direction = solver_>getObjSense(); 

5115  int iColumn; 

5116  const double * dual = cbcRowPrice_; 

5117  down = new double[numberColumns]; 

5118  up = new double[numberColumns]; 

5119  double upSum=1.0e20; 

5120  double downSum = 1.0e20; 

5121  int numberIntegers=0; 

5122  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

5123  CoinBigIndex start = columnStart[iColumn]; 

5124  CoinBigIndex end = start + columnLength[iColumn]; 

5125  double upValue = 0.0; 

5126  double downValue = 0.0; 

5127  double value = direction*objective[iColumn]; 

5128  if (value) { 

5129  if (value>0.0) 

5130  upValue += value; 

5131  else 

5132  downValue = value; 

5133  } 

5134  for (CoinBigIndex j=start;j<end;j++) { 

5135  int iRow = row[j]; 

5136  value = dual[iRow]; 

5137  if (value) { 

5138  value *= element[j]; 

5139  if (value>0.0) 

5140  upValue += value; 

5141  else 

5142  downValue = value; 

5143  } 

5144  } 

5145  // use dj if bigger 

5146  double dj = cbcReducedCost_[iColumn]; 

5147  upValue = CoinMax(upValue,dj); 

5148  downValue = CoinMax(downValue,dj); 

5149  up[iColumn]=upValue; 

5150  down[iColumn]=downValue; 

5151  if (solver_>isInteger(iColumn)) { 

5152  if (!numberNodes_) 

5153  printf("%d  dj %g up %g down %g cost %g\n", 

5154  iColumn,dj,upValue,downValue,objective[iColumn]); 

5155  upSum += upValue; 

5156  downSum += downValue; 

5157  numberIntegers++; 

5158  } 

5159  } 

5160  if (numberIntegers) { 

5161  double smallDown = 0.01*(downSum/((double) numberIntegers)); 

5162  double smallUp = 0.01*(upSum/((double) numberIntegers)); 

5163  for (int i=0;i<numberObjects_;i++) { 

5164  CbcSimpleIntegerDynamicPseudoCost * obj1 = 

5165  dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ; 

5166  if (obj1) { 

5167  iColumn = obj1>columnNumber(); 

5168  double upPseudoCost = obj1>upDynamicPseudoCost(); 

5169  double saveUp = upPseudoCost; 

5170  upPseudoCost = CoinMax(upPseudoCost,smallUp); 

5171  upPseudoCost = CoinMax(upPseudoCost,up[iColumn]); 

5172  upPseudoCost = CoinMax(upPseudoCost,0.1*down[iColumn]); 

5173  obj1>setUpDynamicPseudoCost(upPseudoCost); 

5174  if (upPseudoCost>saveUp&&!numberNodes_) 

5175  printf("For %d up went from %g to %g\n", 

5176  iColumn,saveUp,upPseudoCost); 

5177  double downPseudoCost = obj1>downDynamicPseudoCost(); 

5178  double saveDown = downPseudoCost; 

5179  downPseudoCost = CoinMax(downPseudoCost,smallDown); 

5180  downPseudoCost = CoinMax(downPseudoCost,down[iColumn]); 

5181  downPseudoCost = CoinMax(downPseudoCost,0.1*down[iColumn]); 

5182  obj1>setDownDynamicPseudoCost(downPseudoCost); 

5183  if (downPseudoCost>saveDown&&!numberNodes_) 

5184  printf("For %d down went from %g to %g\n", 

5185  iColumn,saveDown,downPseudoCost); 

5186  } 

5187  } 

5188  } 

5189  delete [] down; 

5190  delete [] up; 

5191  } 

5192  

5193  /* 

5194  Set branching priorities. 

5195  

5196  Setting integer priorities looks pretty robust; the call to findIntegers 

5197  makes sure that SimpleInteger objects are in place. Setting priorities for 

5198  other objects is entirely dependent on their existence, and the routine may 

5199  quietly fail in several directions. 

5200  */ 

5201  

5202  void 

5203  CbcModel::passInPriorities (const int * priorities, 

5204  bool ifObject) 

5205  { 

5206  findIntegers(false); 

5207  int i; 

5208  if (priorities) { 

5209  int i0=0; 

5210  int i1=numberObjects_1; 

5211  if (ifObject) { 

5212  for (i=numberIntegers_;i<numberObjects_;i++) { 

5213  object_[i]>setPriority(priorities[inumberIntegers_]); 

5214  } 

5215  i0=numberIntegers_; 

5216  } else { 

5217  for (i=0;i<numberIntegers_;i++) { 

5218  object_[i]>setPriority(priorities[i]); 

5219  } 

5220  i1=numberIntegers_1; 

5221  } 

5222  messageHandler()>message(CBC_PRIORITY, 

5223  messages()) 

5224  << i0<<i1<<numberObjects_ << CoinMessageEol ; 

5225  } 

5226  } 

5227  

5228  // Delete all object information 

5229  void 

5230  CbcModel::deleteObjects() 

5231  { 

5232  int i; 

5233  for (i=0;i<numberObjects_;i++) 

5234  delete object_[i]; 

5235  delete [] object_; 

5236  object_ = NULL; 

5237  numberObjects_=0; 

5238  findIntegers(true); 

5239  } 

5240  

5241  /*! 

5242  Ensure all attached objects (CbcObjects, heuristics, and cut 

5243  generators) point to this model. 

5244  */ 

5245  void CbcModel::synchronizeModel() 

5246  { 

5247  int i; 

5248  for (i=0;i<numberHeuristics_;i++) 

5249  heuristic_[i]>setModel(this); 

5250  for (i=0;i<numberObjects_;i++) 

5251  object_[i]>setModel(this); 

5252  for (i=0;i<numberCutGenerators_;i++) 

5253  generator_[i]>refreshModel(this); 

5254  } 

5255  

5256  // Fill in integers and create objects 

5257  

5258  /** 

5259  The routine first does a scan to count the number of integer variables. 

5260  It then creates an array, integerVariable_, to store the indices of the 

5261  integer variables, and an array of `objects', one for each variable. 

5262  

5263  The scan is repeated, this time recording the index of each integer 

5264  variable in integerVariable_, and creating an CbcSimpleInteger object that 

5265  contains information about the integer variable. Initially, this is just 

5266  the index and upper & lower bounds. 

5267  

5268  \todo 

5269  Note the assumption in cbc that the first numberIntegers_ objects are 

5270  CbcSimpleInteger. In particular, the code which handles the startAgain 

5271  case assumes that if the object_ array exists it can simply replace the first 

5272  numberInteger_ objects. This is arguably unsafe. 

5273  

5274  I am going to reorder if necessary 

5275  */ 

5276  

5277  void 

5278  CbcModel::findIntegers(bool startAgain,int type) 

5279  { 

5280  assert(solver_); 

5281  /* 

5282  No need to do this if we have previous information, unless forced to start 

5283  over. 

5284  */ 

5285  if (numberIntegers_&&!startAgain&&object_) 

5286  return; 

5287  /* 

5288  Clear out the old integer variable list, then count the number of integer 

5289  variables. 

5290  */ 

5291  delete [] integerVariable_; 

5292  numberIntegers_=0; 

5293  int numberColumns = getNumCols(); 

5294  int iColumn; 

5295  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

5296  if (isInteger(iColumn)) 

5297  numberIntegers_++; 

5298  } 

5299  // Find out how many old noninteger objects there are 

5300  int nObjects=0; 

5301  CbcObject ** oldObject = object_; 

5302  int iObject; 

5303  for (iObject = 0;iObject<numberObjects_;iObject++) { 

5304  CbcSimpleInteger * obj = 

5305  dynamic_cast <CbcSimpleInteger *>(oldObject[iObject]) ; 

5306  if (obj) 

5307  delete oldObject[iObject]; 

5308  else 

5309  oldObject[nObjects++]=oldObject[iObject]; 

5310  } 

5311  

5312  /* 

5313  Found any? Allocate an array to hold the indices of the integer variables. 

5314  Make a large enough array for all objects 

5315  */ 

5316  object_ = new CbcObject * [numberIntegers_+nObjects]; 

5317  numberObjects_=numberIntegers_+nObjects;; 

5318  integerVariable_ = new int [numberIntegers_]; 

5319  /* 

5320  Walk the variables again, filling in the indices and creating objects for 

5321  the integer variables. Initially, the objects hold the index and upper & 

5322  lower bounds. 

5323  */ 

5324  numberIntegers_=0; 

5325  for (iColumn=0;iColumn<numberColumns;iColumn++) { 

5326  if(isInteger(iColumn)) { 

5327  if (!type) 

5328  object_[numberIntegers_] = 

5329  new CbcSimpleInteger(this,numberIntegers_,iColumn); 

5330  else if (type==1) 

5331  object_[numberIntegers_] = 

5332  new CbcSimpleIntegerPseudoCost(this,numberIntegers_,iColumn,0.3); 

5333  integerVariable_[numberIntegers_++]=iColumn; 

5334  } 

5335  } 

5336  // Now append other objects 

5337  memcpy(object_+numberIntegers_,oldObject,nObjects*sizeof(CbcObject *)); 

5338  // Delete old array (just array) 

5339  delete [] oldObject; 

5340  

5341  if (!numberObjects_) 

5342  handler_>message(CBC_NOINT,messages_) << CoinMessageEol ; 

5343  } 

5344  /* If numberBeforeTrust >0 then we are going to use CbcBranchDynamic. 

5345  Scan and convert CbcSimpleInteger objects 

5346  */ 

5347  void 

5348  CbcModel::convertToDynamic() 

5349  { 

5350  int iObject; 

5351  const double * cost = solver_>getObjCoefficients(); 

5352  for (iObject = 0;iObject<numberObjects_;iObject++) { 

5353  CbcSimpleInteger * obj1 = 

5354  dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ; 

5355  CbcSimpleIntegerPseudoCost * obj1a = 

5356  dynamic_cast <CbcSimpleIntegerPseudoCost *>(object_[iObject]) ; 

5357  CbcSimpleIntegerDynamicPseudoCost * obj2 = 

5358  dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[iObject]) ; 

5359  if (obj1&&!obj2) { 

5360  // replace 

5361  int iColumn = obj1>columnNumber(); 

5362  int priority = obj1>priority(); 

5363  int preferredWay = obj1>preferredWay(); 

5364  delete object_[iObject]; 

5365  double costValue = CoinMax(1.0e5,fabs(cost[iColumn])); 

5366  // treat as if will cost what it says up 

5367  double upCost=costValue; 

5368  // and balance at breakeven of 0.3 

5369  double downCost=(0.7*upCost)/0.3; 

5370  if (obj1a) { 

5371  upCost=obj1a>upPseudoCost(); 

5372  downCost=obj1a>downPseudoCost(); 

5373  } 

5374  CbcSimpleIntegerDynamicPseudoCost * newObject = 

5375  new CbcSimpleIntegerDynamicPseudoCost(this,iObject,iColumn,downCost,upCost); 

5376  newObject>setNumberBeforeTrust(numberBeforeTrust_); 

5377  newObject>setPriority(priority); 

5378  newObject>setPreferredWay(preferredWay); 

5379  object_[iObject] = newObject; 

5380  } 

5381  } 

5382  if (branchingMethod_&&(branchingMethod_>whichMethod()&1)==0) { 

5383  // Need a method which can do better 

5384  branchingMethod_=NULL; 

5385  } 

5386  } 

5387  

5388  /* Add in any object information (objects are cloned  owner can delete 

5389  originals */ 

5390  void 

5391  CbcModel::addObjects(int numberObjects, CbcObject ** objects) 

5392  { 

5393  // If integers but not enough objects fudge 

5394  if (numberIntegers_>numberObjects_) 

5395  findIntegers(true); 

5396  /* But if incoming objects inherit from simple integer we just want 

5397  to replace */ 

5398  int numberColumns = solver_>getNumCols(); 

5399  /** mark is 1 if not integer, >=0 if using existing simple integer and 

5400  >=numberColumns if using new integer */ 

5401  int * mark = new int[numberColumns]; 

5402  int i; 

5403  for (i=0;i<numberColumns;i++) 

5404  mark[i]=1; 

5405  int newNumberObjects = numberObjects; 

5406  int newIntegers=0; 

5407  for (i=0;i<numberObjects;i++) { 

5408  CbcSimpleInteger * obj = 

5409  dynamic_cast <CbcSimpleInteger *>(objects[i]) ; 

5410  if (obj) { 

5411  int iColumn = obj>columnNumber(); 

5412  mark[iColumn]=i+numberColumns; 

5413  newIntegers++; 

5414  } 

5415  } 

5416  // and existing 

5417  for (i=0;i<numberObjects_;i++) { 

5418  CbcSimpleInteger * obj = 

5419  dynamic_cast <CbcSimpleInteger *>(object_[i]) ; 

5420  if (obj) { 

5421  int iColumn = obj>columnNumber(); 

5422  if (mark[iColumn]<0) { 

5423  newIntegers++; 

5424  newNumberObjects++; 

5425  mark[iColumn]=i; 

5426  } 

5427  } 

5428  } 

5429  delete [] integerVariable_; 

5430  integerVariable_=NULL; 

5431  if (newIntegers!=numberIntegers_) 

5432  printf("changing number of integers from %d to %d\n", 

5433  numberIntegers_,newIntegers); 

5434  numberIntegers_ = newIntegers; 

5435  integerVariable_ = new int [numberIntegers_]; 

5436  CbcObject ** temp = new CbcObject * [newNumberObjects]; 

5437  // Put integers first 

5438  newIntegers=0; 

5439  numberIntegers_=0; 

5440  for (i=0;i<numberColumns;i++) { 

5441  int which = mark[i]; 

5442  if (which>=0) { 

5443  if (!isInteger(i)) { 

5444  newIntegers++; 

5445  solver_>setInteger(i); 

5446  } 

5447  if (which<numberColumns) { 

5448  temp[numberIntegers_]=object_[which]; 

5449  object_[which]=NULL; 

5450  } else { 

5451  temp[numberIntegers_]=objects[whichnumberColumns]>clone(); 

5452  temp[numberIntegers_]>setModel(this); 

5453  } 

5454  integerVariable_[numberIntegers_++]=i; 

5455  } 

5456  } 

5457  if (newIntegers) 

5458  printf("%d variables were declared integer\n",newIntegers); 

5459  int n=numberIntegers_; 

5460  // Now rest of old 

5461  for (i=0;i<numberObjects_;i++) { 

5462  if (object_[i]) { 

5463  CbcSimpleInteger * obj = 

5464  dynamic_cast <CbcSimpleInteger *>(object_[i]) ; 

5465  if (obj) { 

5466  delete object_[i]; 

5467  } else { 

5468  temp[n++]=object_[i]; 

5469  } 

5470  } 

5471  } 

5472  // and rest of new 

5473  for (i=0;i<numberObjects;i++) { 

5474  CbcSimpleInteger * obj = 

5475  dynamic_cast <CbcSimpleInteger *>(objects[i]) ; 

5476  if (!obj) { 

5477  temp[n]=objects[i]>clone(); 

5478  temp[n++]>setModel(this); 

5479  } 

5480  } 

5481  delete [] mark; 

5482  delete [] object_; 

5483  object_ = temp; 

5484  assert (n==newNumberObjects); 

5485  numberObjects_ = newNumberObjects; 

5486  } 

5487  

5488  /** 

5489  This routine sets the objective cutoff value used for fathoming and 

5490  determining monotonic variables. 

5491  

5492  If the fathoming discipline is strict, a small tolerance is added to the 

5493  new cutoff. This avoids problems due to roundoff when the target value 

5494  is exact. The common example would be an IP with only integer variables in 

5495  the objective. If the target is set to the exact value z of the optimum, 

5496  it's possible to end up fathoming an ancestor of the solution because the 

5497  solver returns z+epsilon. 

5498  

5499  Determining if strict fathoming is needed is best done by analysis. 

5500  In cbc, that's analyseObjective. The default is false. 

5501  

5502  In cbc we always minimize so add epsilon 

5503  */ 

5504  

5505  void CbcModel::setCutoff (double value) 

5506  

5507  { 

5508  #if 0 

5509  double tol = 0 ; 

5510  int fathomStrict = getIntParam(CbcFathomDiscipline) ; 

5511  if (fathomStrict == 1) 

5512  { solver_>getDblParam(OsiDualTolerance,tol) ; 

5513  tol = tol*(1+fabs(value)) ; 

5514  

5515  value += tol ; } 

5516  #endif 

5517  dblParam_[CbcCurrentCutoff]=value; 

5518  // Solvers know about direction 

5519  double direction = solver_>getObjSense(); 

5520  solver_>setDblParam(OsiDualObjectiveLimit,value*direction); } 

5521  

5522  

5523  

5524  /* 

5525  Call this to really test if a valid solution can be feasible. The cutoff is 

5526  passed in as a parameter so that we don't need to worry here after swapping 

5527  solvers. The solution is assumed to be numberColumns in size. If 

5528  fixVariables is true then the bounds of the continuous solver are updated. 

5529  The routine returns the objective value determined by reoptimizing from 

5530  scratch. If the solution is rejected, this will be worse than the cutoff. 

5531  

5532  TODO: There's an issue with getting the correct cutoff value: We update the 

5533  cutoff in the regular solver, but not in continuousSolver_. But our only 

5534  use for continuousSolver_ is verifying candidate solutions. Would it 

5535  make sense to update the cutoff? Then we wouldn't need to step around 

5536  isDualObjectiveLimitReached(). 

5537  */ 

5538  double 

5539  CbcModel::checkSolution (double cutoff, const double *solution, 

5540  bool fixVariables) 

5541  

5542  { 

5543  if (!solverCharacteristics_>solutionAddsCuts()) { 

5544  // Can trust solution 

5545  int numberColumns = solver_>getNumCols(); 

5546  

5547  /* 

5548  Grab the continuous solver (the pristine copy of the problem, made before 

5549  starting to work on the root node). Save the bounds on the variables. 

5550  Install the solution passed as a parameter, and copy it to the model's 

5551  currentSolution_. 

5552  

5553  TODO: This is a beltandsuspenders approach. Once the code has settled 

5554  a bit, we can cast a critical eye here. 

5555  */ 

5556  OsiSolverInterface * saveSolver = solver_; 

5557  if (continuousSolver_) 

5558  solver_ = continuousSolver_; 

5559  // move solution to continuous copy 

5560  solver_>setColSolution(solution); 

5561  // Put current solution in safe place 

5562  // Point to current solution 

5563  const double * save = testSolution_; 

5564  // Safe as will be const inside infeasibility() 

5565  testSolution_ = solver_>getColSolution(); 

5566  //memcpy(currentSolution_,solver_>getColSolution(), 

5567  // numberColumns*sizeof(double)); 

5568  //solver_>messageHandler()>setLogLevel(4); 

5569  

5570  // save original bounds 

5571  double * saveUpper = new double[numberColumns]; 

5572  double * saveLower = new double[numberColumns]; 

5573  memcpy(saveUpper,getColUpper(),numberColumns*sizeof(double)); 

5574  memcpy(saveLower,getColLower(),numberColumns*sizeof(double)); 

5575  

5576  /* 

5577  Run through the objects and use feasibleRegion() to set variable bounds 

5578  so as to fix the variables specified in the objects at their value in this 

5579  solution. Since the object list contains (at least) one object for every 

5580  integer variable, this has the effect of fixing all integer variables. 

5581  */ 

5582  int i; 

5583  for (i=0;i<numberObjects_;i++) 

5584  object_[i]>feasibleRegion(); 

5585  // We can switch off check 

5586  if ((specialOptions_&4)==0) { 

5587  if ((specialOptions_&2)==0&&solverCharacteristics_>warmStart()) { 

5588  /* 

5589  Remove any existing warm start information to be sure there is no 

5590  residual influence on initialSolve(). 

5591  */ 

5592  CoinWarmStartBasis *slack = 

5593  dynamic_cast<CoinWarmStartBasis *>(solver_>getEmptyWarmStart()) ; 

5594  solver_>setWarmStart(slack); 

5595  delete slack ; 

5596  } 

5597  // Give a hint to do dual 

5598  bool saveTakeHint; 

5599  OsiHintStrength saveStrength; 

5600  bool gotHint = (solver_>getHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength)); 

5601  assert (gotHint); 

5602  solver_>setHintParam(OsiDoDualInInitial,true,OsiHintTry); 

5603  solver_>initialSolve(); 

5604  if (!solver_>isProvenOptimal()) 

5605  { //printf("checkSolution infeas! Retrying wihout scaling.\n"); 

5606  bool saveTakeHint; 

5607  OsiHintStrength saveStrength; 

5608  bool savePrintHint; 

5609  //solver_>writeMps("infeas"); 

5610  bool gotHint = (solver_>getHintParam(OsiDoReducePrint,savePrintHint,saveStrength)); 

5611  gotHint = (solver_>getHintParam(OsiDoScale,saveTakeHint,saveStrength)); 

5612  solver_>setHintParam(OsiDoScale,false,OsiHintTry); 

5613  //solver_>setHintParam(OsiDoReducePrint,false,OsiHintTry) ; 

5614  solver_>setHintParam(OsiDoDualInInitial,false,OsiHintTry); 

5615  solver_>initialSolve(); 

5616  solver_>setHintParam(OsiDoScale,saveTakeHint,saveStrength); 

5617  //solver_>setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ; 

5618  } 

5619  //assert(solver_>isProvenOptimal()); 

5620  solver_>setHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength); 

5621  } 

5622  double objectiveValue = solver_>getObjValue()*solver_>getObjSense(); 

5623  

5624  /* 

5625  Check that the solution still beats the objective cutoff. 

5626  

5627  If it passes, make a copy of the primal variable values and do some 

5628  cleanup and checks: 

5629  + Values of all variables are are within original bounds and values of 

5630  all integer variables are within tolerance of integral. 

5631  + There are no constraint violations. 

5632  There really should be no need for the check against original bounds. 

5633  Perhaps an opportunity for a sanity check? 

5634  */ 

5635  if ((solver_>isProvenOptimal()(specialOptions_&4)!=0) && objectiveValue <= cutoff) { 

5636  double * solution = new double[numberColumns]; 

5637  memcpy(solution ,solver_>getColSolution(),numberColumns*sizeof(double)) ; 

5638  

5639  int iColumn; 

5640  #ifndef NDEBUG 

5641  double integerTolerance = getIntegerTolerance() ; 

5642  #endif 

5643  for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) { 

5644  double value = solution[iColumn] ; 

5645  value = CoinMax(value, saveLower[iColumn]) ; 

5646  value = CoinMin(value, saveUpper[iColumn]) ; 

5647  if (solver_>isInteger(iColumn)) 

5648  assert(fabs(valuesolution[iColumn]) <= integerTolerance) ; 

5649  solution[iColumn] = value ; 

5650  } 

5651  if ((specialOptions_&16)==0) { 

5652  const double * rowLower = solver_>getRowLower() ; 

5653  const double * rowUpper = solver_>getRowUpper() ; 

5654  int numberRows = solver_>getNumRows() ; 

5655  double *rowActivity = new double[numberRows] ; 

5656  memset(rowActivity,0,numberRows*sizeof(double)) ; 

5657  solver_>getMatrixByCol()>times(solution,rowActivity) ; 

5658  double primalTolerance ; 

5659  solver_>getDblParam(OsiPrimalTolerance,primalTolerance) ; 

5660  double largestInfeasibility =0.0; 

5661  for (i=0 ; i < numberRows ; i++) { 

5662  largestInfeasibility = CoinMax(largestInfeasibility, 

5663  rowLower[i]rowActivity[i]); 

5664  largestInfeasibility = CoinMax(largestInfeasibility, 

5665  rowActivity[i]rowUpper[i]); 

5666  } 

5667  if (largestInfeasibility>100.0*primalTolerance) { 

5668  handler_>message(CBC_NOTFEAS3, messages_) 

5669  << largestInfeasibility << CoinMessageEol ; 

5670  objectiveValue=1.0e50 ; 

5671  } 

5672  delete [] rowActivity ; 

5673  } 

5674  delete [] solution; 

5675  } else { 

5676  objectiveValue=1.0e50 ; 

5677  } 

5678  /* 

5679  Regardless of what we think of the solution, we may need to restore the 

5680  original bounds of the continuous solver. Unfortunately, const'ness 

5681  prevents us from simply reversing the memcpy used to make these snapshots. 

5682  */ 

5683  if (!fixVariables) 

5684  { for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) 

5685  { solver_>setColLower(iColumn,saveLower[iColumn]) ; 

5686  solver_>setColUpper(iColumn,saveUpper[iColumn]) ; } } 

5687  delete [] saveLower; 

5688  delete [] saveUpper; 

5689  

5690  /* 

5691  Restore the usual solver. 

5692  */ 

5693  solver_=saveSolver; 

5694  testSolution_ = save; 

5695  return objectiveValue; 

5696  } else { 

5697  // Outer approximation or similar 

5698  //If this is true then the solution comes from the nlp we don't need to resolve the same nlp with ipopt 

5699  //solverCharacteristics_>setSolver(solver_); 

5700  bool solutionComesFromNlp = solverCharacteristics_>bestObjectiveValue()<cutoff; 

5701  double objectiveValue; 

5702  int numberColumns = solver_>getNumCols(); 

5703  double *saveLower = NULL; 

5704  double * saveUpper = NULL; 

5705  

5706  if(! solutionComesFromNlp)//Otherwise solution already comes from ipopt and cuts are known 

5707  { 

5708  if(fixVariables)//Will temporarily fix all integer valued var 

5709  { 

5710  // save original bounds 

5711  saveUpper = new double[numberColumns]; 

5712  saveLower = new double[numberColumns]; 

5713  memcpy(saveUpper,solver_>getColUpper(),numberColumns*sizeof(double)); 

5714  memcpy(saveLower,solver_>getColLower(),numberColumns*sizeof(double)); 

5715  //in any case solution should be already loaded into solver_ 

5716  /* 

5717  Run through the objects and use feasibleRegion() to set variable bounds 

5718  so as to fix the variables specified in the objects at their value in this 

5719  solution. Since the object list contains (at least) one object for every 

5720  integer variable, this has the effect of fixing all integer variables. 

5721  */ 

5722  const double * save = testSolution_; 

5723  testSolution_ = solution; 

5724  for (int i=0;i<numberObjects_;i++) 

5725  object_[i]>feasibleRegion(); 

5726  testSolution_ = save; 

5727  solver_>resolve(); 

5728  } 

5729  

5730  /* 

5731  Now step through the cut generators and see if any of them are flagged to 

5732  run when a new solution is discovered. Only global cuts are useful. 

5733  (The solution being evaluated may not correspond to the current location in the 

5734  search tree  discovered by heuristic, for example.) 

5735  */ 

5736  OsiCuts theseCuts; 

5737  int i; 

5738  int lastNumberCuts=0; 

5739  for (i=0;i<numberCutGenerators_;i++) 

5740  { 

5741  if (generator_[i]>atSolution()) 

5742  { 

5743  generator_[i]>generateCuts(theseCuts,true,NULL); 

5744  int numberCuts = theseCuts.sizeRowCuts(); 

5745  for (int j=lastNumberCuts;j<numberCuts;j++) 

5746  { 

5747  const OsiRowCut * thisCut = theseCuts.rowCutPtr(j); 

5748  if (thisCut>globallyValid()) 

5749  { 

5750  // if ((specialOptions_&1)!=0) 

5751  // { 

5752  // /* As these are global cuts  

5753  // a) Always get debugger object 

5754  // b) Not fatal error to cutoff optimal (if we have just got optimal) 

5755  // */ 

5756  // const OsiRowCutDebugger *debugger = solver_>getRowCutDebuggerAlways() ; 

5757  // if (debugger) 

5758  // { 

5759  // if(debugger>invalidCut(*thisCut)) 

5760  // printf("ZZZZ Global cut  cuts off optimal solution!\n"); 

5761  // } 

5762  // } 

5763  // add to global list 

5764  globalCuts_.insert(*thisCut); 

5765  generator_[i]>incrementNumberCutsInTotal(); 

5766  } 

5767  } 

5768  } 

5769  } 

5770  // int numberCuts = theseCuts.sizeColCuts(); 

5771  // for (i=0;i<numberCuts;i++) { 

5772  // const OsiColCut * thisCut = theseCuts.colCutPtr(i); 

5773  // if (thisCut>globallyValid()) { 

5774  // // add to global list 

5775  // globalCuts_.insert(*thisCut); 

5776  // } 

5777  // } 

5778  //have to retrieve the solution and its value from the nlp 

5779  } 

5780  double newObjectiveValue=cutoff; 

5781  if(solverCharacteristics_>solution(newObjectiveValue, 

5782  const_cast<double *> (solution), 

5783  numberColumns)) 

5784  { 

5785  objectiveValue = newObjectiveValue; 

5786  } 

5787  else 

5788  { 

5789  objectiveValue = 2e50; 

5790  } 

5791  if (!solutionComesFromNlp && fixVariables) 

5792  { 

5793  for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) 

5794  { 

5795  solver_>setColLower(iColumn,saveLower[iColumn]) ; 

5796  solver_>setColUpper(iColumn,saveUpper[iColumn]) ; 

5797  } 

5798  delete [] saveLower; 

5799  delete [] saveUpper; 

5800  } 

5801  //If the variables were fixed the cutting plane procedure may have believed that the node could be fathomed 

5802  //reestablish truth. should do no harm for non nlp 

5803  if(!solutionComesFromNlp && fixVariables) 

5804  solverCharacteristics_>setMipBound(COIN_DBL_MAX); 

5805  return objectiveValue; 

5806  } 

5807  } 

5808  

5809  /* 

5810  Call this routine from anywhere when a solution is found. The solution 

5811  vector is assumed to contain one value for each structural variable. 

5812  

5813  The first action is to run checkSolution() to confirm the objective and 

5814  feasibility. If this check causes the solution to be rejected, we're done. 

5815  If fixVariables = true, the variable bounds held by the continuous solver 

5816  will be left fixed to the values in the solution; otherwise they are 

5817  restored to the original values. 

5818  

5819  If the solution is accepted, install it as the best solution. 

5820  

5821  The routine also contains a hook to run any cut generators that are flagged 

5822  to run when a new solution is discovered. There's a potential hazard because 

5823  the cut generators see the continuous solver >after< possible restoration of 

5824  original bounds (which may well invalidate the solution). 

5825  */ 

5826  

5827  void 

5828  CbcModel::setBestSolution (CBC_Message how, 

5829  double & objectiveValue, const double *solution, 

5830  bool fixVariables) 

5831  

5832  { 

5833  if (!solverCharacteristics_>solutionAddsCuts()) { 

5834  // Can trust solution 

5835  double cutoff = CoinMin(getCutoff(),bestObjective_) ; 

5836  

5837  /* 

5838  Double check the solution to catch pretenders. 

5839  */ 

5840  objectiveValue = checkSolution(cutoff,solution,fixVariables); 

5841  if (objectiveValue > cutoff) { 

5842  if (objectiveValue>1.0e30) 

5843  handler_>message(CBC_NOTFEAS1, messages_) << CoinMessageEol ; 

5844  else 

5845  handler_>message(CBC_NOTFEAS2, messages_) 

5846  << objectiveValue << cutoff << CoinMessageEol ; 

5847  } else { 


