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  #include <cassert> 

12  #include <cmath> 

13  #include <cfloat> 

14  

15  #include "OsiSolverInterface.hpp" 

16  #include "CoinWarmStartBasis.hpp" 

17  #include "CoinPackedMatrix.hpp" 

18  #include "CoinHelperFunctions.hpp" 

19  #include "CbcBranchActual.hpp" 

20  #include "CbcBranchDynamic.hpp" 

21  #include "CbcHeuristic.hpp" 

22  #include "CbcModel.hpp" 

23  #include "CbcStatistics.hpp" 

24  #include "CbcStrategy.hpp" 

25  #include "CbcMessage.hpp" 

26  #include "OsiRowCut.hpp" 

27  #include "OsiColCut.hpp" 

28  #include "OsiRowCutDebugger.hpp" 

29  #include "OsiCuts.hpp" 

30  #include "CbcCountRowCut.hpp" 

31  #include "CbcCutGenerator.hpp" 

32  // include Probing 

33  #include "CglProbing.hpp" 

34  

35  #define COIN_USE_CLP 

36  #ifdef COIN_USE_CLP 

37  // include Presolve from Clp 

38  #include "ClpPresolve.hpp" 

39  #include "OsiClpSolverInterface.hpp" 

40  #include "ClpEventHandler.hpp" 

41  #endif 

42  

43  #include "CoinTime.hpp" 

44  

45  #include "CbcCompareActual.hpp" 

46  #include "CbcTree.hpp" 

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

48  

49  namespace { 

50  

51  // 

52  // Returns the greatest common denominator of two 

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

54  // 

55  static int gcd(int a, int b) 

56  { 

57  int remainder = 1; 

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

59  if(a > b) { 

60  // Swap a and b 

61  int temp = a; 

62  a = b; 

63  b = temp; 

64  } 

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

66  if (!a) { 

67  if (b) { 

68  return b; 

69  } else { 

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

71  abort(); 

72  } 

73  } 

74  while (remainder) { 

75  remainder = b % a; 

76  b = a; 

77  a = remainder; 

78  } 

79  return b; 

80  } 

81  

82  

83  

84  #ifdef CHECK_NODE_FULL 

85  

86  /* 

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

88  is 

89  

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

91  

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

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

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

95  made to see if the numberPointingToThis_ field agrees. 

96  */ 

97  

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

99  

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

101  

102  int j ; 

103  int nNodes = branchingTree>size() ; 

104  # define MAXINFO 1000 

105  int *count = new int [MAXINFO] ; 

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

107  int nInfo = 0 ; 

108  /* 

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

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

111  branches remaining. 

112  

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

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

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

116  common ancestor. 

117  */ 

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

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

120  CbcNodeInfo *nodeInfo = node>nodeInfo() ; 

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

122  assert(change) ; 

123  while (nodeInfo) 

124  { int k ; 

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

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

127  if (k == nInfo) 

128  { assert(nInfo < MAXINFO) ; 

129  nInfo++ ; 

130  info[k] = nodeInfo ; 

131  count[k] = 0 ; } 

132  nodeInfo = nodeInfo>parent() ; } } 

133  /* 

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

135  increment the corresponding count. 

136  */ 

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

138  { CbcNodeInfo *nodeInfo = info[j] ; 

139  nodeInfo = nodeInfo>parent() ; 

140  if (nodeInfo) 

141  { int k ; 

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

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

144  assert (k < nInfo) ; 

145  count[k]++ ; } } 

146  /* 

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

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

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

150  references (unexplored branches, numberBranchesLeft()). 

151  */ 

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

153  CbcNodeInfo * nodeInfo = info[j] ; 

154  if (nodeInfo) { 

155  int k ; 

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

157  if (nodeInfo == info[k]) 

158  break ; 

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

160  nodeInfo, 

161  nodeInfo>numberBranchesLeft(), 

162  nodeInfo>numberPointingToThis()) ; 

163  assert(nodeInfo>numberPointingToThis() == 

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

165  

166  delete [] count ; 

167  delete [] info ; 

168  

169  return ; } 

170  

171  #endif /* CHECK_NODE_FULL */ 

172  

173  

174  

175  #ifdef CHECK_CUT_COUNTS 

176  

177  /* 

178  Routine to verify that cut reference counts are correct. 

179  */ 

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

181  

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

183  

184  int j ; 

185  int nNodes = branchingTree>size() ; 

186  

187  /* 

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

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

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

191  */ 

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

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

194  CbcNodeInfo * nodeInfo = node>nodeInfo() ; 

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

196  while (nodeInfo) 

197  { int k ; 

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

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

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

201  nodeInfo = nodeInfo>parent() ; } } 

202  /* 

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

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

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

206  */ 

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

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

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

210  CbcNodeInfo *nodeInfo = node>nodeInfo(); 

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

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

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

214  node>objectiveValue()) ; 

215  

216  model.addCuts1(node,debugws) ; 

217  

218  int i ; 

219  int numberRowsAtContinuous = model.numberRowsAtContinuous() ; 

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

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

222  { CoinWarmStartBasis::Status status = 

223  debugws>getArtifStatus(i+numberRowsAtContinuous) ; 

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

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

226  

227  while (nodeInfo) 

228  { nodeInfo = nodeInfo>parent() ; 

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

230  printf("\n") ; 

231  delete debugws ; } 

232  /* 

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

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

235  

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

237  */ 

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

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

240  CbcNodeInfo *nodeInfo = node>nodeInfo(); 

241  while (nodeInfo) 

242  { int k ; 

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

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

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

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

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

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

249  else 

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

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

252  cut>tempNumber_ = 1 ; } } 

253  nodeInfo = nodeInfo>parent() ; } } 

254  

255  return ; } 

256  

257  #endif /* CHECK_CUT_COUNTS */ 

258  

259  } 

260  

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

262  

263  

264  

265  void 

266  CbcModel::analyzeObjective () 

267  /* 

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

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

270  and grabs the largest objective coefficient associated with an unfixed 

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

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

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

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

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

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

277  

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

279  

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

281  */ 

282  { const double *objective = getObjCoefficients() ; 

283  const double *lower = getColLower() ; 

284  const double *upper = getColUpper() ; 

285  /* 

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

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

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

289  

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

291  fathoming discipline to strict. 

292  */ 

293  double maximumCost = 0.0 ; 

294  bool possibleMultiple = true ; 

295  int iColumn ; 

296  int numberColumns = getNumCols() ; 

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

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

299  { if (isInteger(iColumn)) 

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

301  else if (objective[iColumn]) 

302  possibleMultiple = false ; } } 

303  setIntParam(CbcModel::CbcFathomDiscipline,possibleMultiple) ; 

304  /* 

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

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

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

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

309  later). 

310  

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

312  */ 

313  if (possibleMultiple&&maximumCost) 

314  { int increment = 0 ; 

315  double multiplier = 2520.0 ; 

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

317  multiplier *= 10.0 ; 

318  

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

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

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

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

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

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

325  { increment = 0 ; 

326  break ; } 

327  else if (!increment) 

328  { increment = nearest ; } 

329  else 

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

331  /* 

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

333  */ 

334  if (increment) 

335  { double value = increment ; 

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

337  value /= multiplier ; 

338  if (value*0.999 > cutoff) 

339  { messageHandler()>message(CBC_INTEGERINCREMENT, 

340  messages()) 

341  << value << CoinMessageEol ; 

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

343  

344  return ; 

345  } 

346  

347  

348  

349  /** 

350  \todo 

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

352  routine, after calling the solver for an initial solution 

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

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

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

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

357  checked. 

358  

359  */ 

360  

361  /* 

362  The overall flow can be divided into three stages: 

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

364  do all the setup for B&C. 

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

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

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

368  

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

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

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

372  system held by the solver. 

373  */ 

374  

375  void CbcModel::branchAndBound(int doStatistics) 

376  

377  { 

378  // Set up strategies 

379  if (strategy_) { 

380  strategy_>setupCutGenerators(*this); 

381  strategy_>setupHeuristics(*this); 

382  // Set strategy print level to models 

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

384  strategy_>setupOther(*this); 

385  } 

386  bool eventHappened=false; 

387  ClpEventHandler * eventHandler=NULL; 

388  #ifdef COIN_USE_CLP 

389  { 

390  OsiClpSolverInterface * clpSolver 

391  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

392  if (clpSolver) { 

393  ClpSimplex * clpSimplex = clpSolver>getModelPtr(); 

394  eventHandler = clpSimplex>eventHandler(); 

395  } 

396  } 

397  #endif 

398  if (!nodeCompare_) 

399  nodeCompare_=new CbcCompareDefault();; 

400  # ifdef CBC_DEBUG 

401  std::string problemName ; 

402  solver_>getStrParam(OsiProbName,problemName) ; 

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

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

405  # endif 

406  /* 

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

408  */ 

409  status_ = 0 ; 

410  secondaryStatus_ = 0; 

411  phase_=0; 

412  /* 

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

414  CbcSimpleInteger object for each integer variable. 

415  */ 

416  findIntegers(false) ; 

417  // If dynamic pseudo costs then do 

418  if (numberBeforeTrust_>0) 

419  convertToDynamic(); 

420  

421  /* 

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

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

424  */ 

425  synchronizeModel() ; 

426  /* 

427  Capture a time stamp before we start. 

428  */ 

429  dblParam_[CbcStartSeconds] = CoinCpuTime(); 

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

431  continuousObjective_ = COIN_DBL_MAX ; 

432  /* 

433  Solve the relaxation. 

434  

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

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

437  continuous relaxation. 

438  */ 

439  bool feasible = resolve() ; 

440  /* 

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

442  continue with processing the root node. 

443  */ 

444  if (!feasible) 

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

446  status_ = 0 ; 

447  secondaryStatus_ = 1; 

448  originalContinuousObjective_ = COIN_DBL_MAX; 

449  return ; } 

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

451  originalContinuousObjective_ = solver_>getObjValue(); 

452  bestPossibleObjective_=originalContinuousObjective_; 

453  sumChangeObjective1_=0.0; 

454  sumChangeObjective2_=0.0; 

455  /* 

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

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

458  see if it cuts off the optimal answer. 

459  */ 

460  // If debugger exists set specialOptions_ bit 

461  if (solver_>getRowCutDebuggerAlways()) 

462  specialOptions_ = 1; 

463  

464  # ifdef CBC_DEBUG 

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

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

467  if (solver_>getRowCutDebuggerAlways()) 

468  specialOptions_ = 1; 

469  # endif 

470  

471  /* 

472  Begin setup to process a feasible root node. 

473  */ 

474  bestObjective_ = CoinMin(bestObjective_,1.0e50) ; 

475  numberSolutions_ = 0 ; 

476  stateOfSearch_=0; 

477  numberHeuristicSolutions_ = 0 ; 

478  // Everything is minimization 

479  double cutoff=getCutoff() ; 

480  double direction = solver_>getObjSense() ; 

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

482  messageHandler()>message(CBC_CUTOFF_WARNING1, 

483  messages()) 

484  << cutoff << cutoff << CoinMessageEol ; 

485  if (cutoff > bestObjective_) 

486  cutoff = bestObjective_ ; 

487  setCutoff(cutoff) ; 

488  /* 

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

490  */ 

491  int numberColumns = getNumCols() ; 

492  if (!currentSolution_) 

493  currentSolution_ = new double[numberColumns] ; 

494  testSolution_ = currentSolution_; 

495  /* 

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

497  constraint system (aka the continuous system). 

498  */ 

499  continuousSolver_ = solver_>clone() ; 

500  #ifdef COIN_USE_CLP 

501  { 

502  OsiClpSolverInterface * clpSolver 

503  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

504  if (clpSolver) { 

505  ClpSimplex * clpSimplex = clpSolver>getModelPtr(); 

506  // take off names 

507  clpSimplex>dropNames(); 

508  } 

509  } 

510  #endif 

511  

512  numberRowsAtContinuous_ = getNumRows() ; 

513  /* 

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

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

516  installed. 

517  */ 

518  analyzeObjective() ; 

519  /* 

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

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

522  that produced a given cut. 

523  */ 

524  int maximumWhich = 1000 ; 

525  int * whichGenerator = new int[maximumWhich] ; 

526  int currentNumberCuts = 0 ; 

527  maximumNumberCuts_ = 0 ; 

528  currentNumberCuts_ = 0 ; 

529  delete [] addedCuts_ ; 

530  addedCuts_ = NULL ; 

531  /* 

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

533  (problems which require further exploration). 

534  */ 

535  tree_>setComparison(*nodeCompare_) ; 

536  /* 

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

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

539  */ 

540  maximumDepth_ = 10 ; 

541  delete [] walkback_ ; 

542  walkback_ = new CbcNodeInfo * [maximumDepth_] ; 

543  /* 

544  Used to generate bound edits for CbcPartialNodeInfo. 

545  */ 

546  double * lowerBefore = new double [numberColumns] ; 

547  double * upperBefore = new double [numberColumns] ; 

548  /* 

549  

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

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

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

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

554  

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

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

557  removed any slack cuts (numberOldActiveCuts and numberNewCuts have been 

558  adjusted accordingly). 

559  

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

561  

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

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

564  */ 

565  phase_=1; 

566  int iCutGenerator; 

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

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

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

570  } 

571  OsiCuts cuts ; 

572  int anyAction = 1 ; 

573  int numberOldActiveCuts = 0 ; 

574  int numberNewCuts = 0 ; 

575  // Array to mark solution 

576  delete [] usedInSolution_; 

577  usedInSolution_ = new int[numberColumns]; 

578  CoinZeroN(usedInSolution_,numberColumns); 

579  /* 

580  For printing totals and for CbcNode (numberNodes_) 

581  */ 

582  numberIterations_ = 0 ; 

583  numberNodes_ = 0 ; 

584  numberNodes2_ = 0 ; 

585  int maximumStatistics=0; 

586  CbcStatistics ** statistics = NULL; 

587  // Do on switch 

588  if (doStatistics) { 

589  maximumStatistics=10000; 

590  statistics = new CbcStatistics * [maximumStatistics]; 

591  memset(statistics,0,maximumStatistics*sizeof(CbcStatistics *)); 

592  } 

593  

594  { int iObject ; 

595  int preferredWay ; 

596  int numberUnsatisfied = 0 ; 

597  memcpy(currentSolution_,solver_>getColSolution(), 

598  numberColumns*sizeof(double)) ; 

599  

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

601  { double infeasibility = 

602  object_[iObject]>infeasibility(preferredWay) ; 

603  if (infeasibility ) numberUnsatisfied++ ; } 

604  if (numberUnsatisfied) { 

605  feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_, 

606  NULL,numberOldActiveCuts,numberNewCuts, 

607  maximumWhich, whichGenerator) ; 

608  } 

609  } 

610  // make cut generators less aggressive 

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

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

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

614  } 

615  currentNumberCuts_ = numberNewCuts ; 

616  /* 

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

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

619  currently uses strong branching to evaluate branch object candidates, 

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

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

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

623  

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

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

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

627  

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

629  */ 

630  

631  bool resolved = false ; 

632  CbcNode *newNode = NULL ; 

633  if (feasible) 

634  { newNode = new CbcNode ; 

635  newNode>setObjectiveValue(direction*solver_>getObjValue()) ; 

636  anyAction = 1 ; 

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

638  CbcNode fakeNode; 

639  if (!currentNode_) { 

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

641  currentNode_=&fakeNode; 

642  } 

643  phase_=3; 

644  // only allow twenty passes 

645  int numberPassesLeft=20; 

646  while (anyAction == 1) 

647  { 

648  if (numberBeforeTrust_<=0 ) { 

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

650  } else { 

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

652  if (anyAction==3) 

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

654  } 

655  numberPassesLeft; 

656  if (anyAction == 1) 

657  { feasible = resolve() ; 

658  resolved = true ; 

659  # ifdef CBC_DEBUG 

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

661  solver_>getObjValue(), 

662  solver_>getNumRows()) ; 

663  # endif 

664  if (!feasible) anyAction = 2 ; } 

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

666  { delete newNode ; 

667  newNode = NULL ; 

668  feasible = false ; } } } 

669  /* 

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

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

672  current objective cutoff. 

673  */ 

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

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

676  CbcNode * rootNode = newNode; 

677  /* 

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

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

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

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

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

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

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

685  a valid solution for use by setBestSolution(). 

686  */ 

687  CoinWarmStartBasis *lastws = 0 ; 

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

689  { if (resolved) 

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

691  takeOffCuts(cuts,whichGenerator,numberOldActiveCuts,numberNewCuts, 

692  needValidSolution) ; 

693  # ifdef CHECK_CUT_COUNTS 

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

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

696  numberRowsAtContinuous_+numberNewCuts+numberOldActiveCuts) ; 

697  const CoinWarmStartBasis* debugws = 

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

699  debugws>print() ; 

700  delete debugws ; } 

701  # endif 

702  } 

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

704  newNode>nodeInfo()>addCuts(cuts, 

705  newNode>numberBranches(),whichGenerator) ; 

706  /* 

707  Courtesy of createInfo, there's now a full basis stashed in 

708  newNode>nodeInfo_>basis_. We're about to make two more copies, lastws and 

709  model.basis_. 

710  

711  (jf) With some thought I should be able to get rid of lastws and use 

712  basis_. 

713  (lh) I agree, but haven't pursued it to the end. 

714  */ 

715  if (basis_) delete basis_ ; 

716  basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

717  if (lastws) delete lastws ; 

718  lastws = dynamic_cast<CoinWarmStartBasis*>(basis_>clone()) ; } 

719  /* 

720  Continuous data to be used later 

721  */ 

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

723  continuousInfeasibilities_ = 0 ; 

724  if (newNode) 

725  { continuousObjective_ = newNode>objectiveValue() ; 

726  delete [] continuousSolution_; 

727  continuousSolution_ = CoinCopyOfArray(solver_>getColSolution(), 

728  numberColumns); 

729  continuousInfeasibilities_ = newNode>numberUnsatisfied() ; } 

730  /* 

731  Bound may have changed so reset in objects 

732  */ 

733  { int i ; 

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

735  object_[i]>resetBounds() ; } 

736  bool stoppedOnGap = false ; 

737  /* 

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

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

740  integral solution at the root node. 

741  

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

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

744  */ 

745  double bestValue = 0.0 ; 

746  if (newNode) { 

747  bestValue = newNode>objectiveValue(); 

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

749  newNode>initializeInfo() ; 

750  tree_>push(newNode) ; 

751  if (statistics) { 

752  if (numberNodes2_==maximumStatistics) { 

753  maximumStatistics = 2*maximumStatistics; 

754  CbcStatistics ** temp = new CbcStatistics * [maximumStatistics]; 

755  memset(temp,0,maximumStatistics*sizeof(CbcStatistics *)); 

756  memcpy(temp,statistics,numberNodes2_*sizeof(CbcStatistics *)); 

757  delete [] statistics; 

758  statistics=temp; 

759  } 

760  assert (!statistics[numberNodes2_]); 

761  statistics[numberNodes2_]=new CbcStatistics(newNode); 

762  } 

763  numberNodes2_++; 

764  # ifdef CHECK_NODE 

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

766  # endif 

767  } else { 

768  // continuous is integer 

769  double objectiveValue = newNode>objectiveValue(); 

770  setBestSolution(CBC_SOLUTION,objectiveValue, 

771  solver_>getColSolution()) ; 

772  delete newNode ; 

773  newNode = NULL ; 

774  } 

775  } 

776  

777  if (printFrequency_ <= 0) { 

778  printFrequency_ = 1000 ; 

779  if (getNumCols() > 2000) 

780  printFrequency_ = 100 ; 

781  } 

782  /* Tell solver we are in Branch and Cut 

783  Could use last parameter for subtle differences */ 

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

785  /* 

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

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

788  */ 

789  int numberLongStrong=0; 

790  /* 

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

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

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

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

795  add the node to the live set. 

796  

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

798  than the current objective cutoff. 

799  */ 

800  while (!tree_>empty()) 

801  { if (cutoff > getCutoff()) { 

802  if (eventHandler) { 

803  if (!eventHandler>event(ClpEventHandler::solution)) { 

804  eventHappened=true; // exit 

805  } 

806  } 

807  // Do from deepest 

808  tree_>cleanTree(this, getCutoff(),bestPossibleObjective_) ; 

809  nodeCompare_>newSolution(this) ; 

810  nodeCompare_>newSolution(this,continuousObjective_, 

811  continuousInfeasibilities_) ; 

812  tree_>setComparison(*nodeCompare_) ; 

813  if (tree_>empty()) 

814  break; // finished 

815  } 

816  cutoff = getCutoff() ; 

817  /* 

818  Periodic activities: Opportunities to 

819  + tweak the nodeCompare criteria, 

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

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

822  */ 

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

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

825  // redo tree if wanted 

826  if (redoTree) 

827  tree_>setComparison(*nodeCompare_) ; 

828  } 

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

830  int j ; 

831  int nNodes = tree_>size() ; 

832  bestPossibleObjective_ = 1.0e100 ; 

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

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

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

836  bestPossibleObjective_ = node>objectiveValue() ; 

837  } 

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

839  << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_ 

840  << CoinMessageEol ; 

841  if (eventHandler) { 

842  if (!eventHandler>event(ClpEventHandler::treeStatus)) { 

843  eventHappened=true; // exit 

844  } 

845  } 

846  } 

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

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

849  stateOfSearch_=3; 

850  // See if can stop on gap 

851  double testGap = CoinMax(dblParam_[CbcAllowableGap], 

852  CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_)) 

853  *dblParam_[CbcAllowableFractionGap]); 

854  if (bestObjective_bestPossibleObjective_ < testGap) { 

855  stoppedOnGap = true ; 

856  } 

857  

858  # ifdef CHECK_NODE_FULL 

859  verifyTreeNodes(tree_,*this) ; 

860  # endif 

861  # ifdef CHECK_CUT_COUNTS 

862  verifyCutCounts(tree_,*this) ; 

863  # endif 

864  

865  /* 

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

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

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

869  active subproblem. 

870  */ 

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

872  // Possible one on tree worse than cutoff 

873  if (!node) 

874  continue; 

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

876  #ifdef CBC_DEBUG 

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

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

879  node>guessedObjectiveValue()); 

880  #endif 

881  // Save clone in branching decision 

882  if(branchingMethod_) 

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

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

885  // Say not on optimal path 

886  bool onOptimalPath=false; 

887  # ifdef CHECK_NODE 

888  /* 

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

890  branching object is something other than an integer variable. 

891  This needs some thought. 

892  */ 

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

894  node>nodeInfo()>numberBranchesLeft(), 

895  node>nodeInfo()>numberPointingToThis()) ; 

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

897  node>depth(),node>objectiveValue(), 

898  node>numberUnsatisfied(), 

899  integerVariable_[node>variable()]) ; 

900  # endif 

901  

902  /* 

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

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

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

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

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

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

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

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

911  live set. 

912  */ 

913  CbcNodeInfo * nodeInfo = node>nodeInfo() ; 

914  newNode = NULL ; 

915  if (!addCuts(node,lastws)) 

916  { int i ; 

917  const double * lower = getColLower() ; 

918  const double * upper = getColUpper() ; 

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

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

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

922  bool deleteNode ; 

923  if (node>branch()) 

924  { 

925  // set nodenumber correctly 

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

927  tree_>push(node) ; 

928  if (statistics) { 

929  if (numberNodes2_==maximumStatistics) { 

930  maximumStatistics = 2*maximumStatistics; 

931  CbcStatistics ** temp = new CbcStatistics * [maximumStatistics]; 

932  memset(temp,0,maximumStatistics*sizeof(CbcStatistics *)); 

933  memcpy(temp,statistics,numberNodes2_*sizeof(CbcStatistics *)); 

934  delete [] statistics; 

935  statistics=temp; 

936  } 

937  assert (!statistics[numberNodes2_]); 

938  statistics[numberNodes2_]=new CbcStatistics(node); 

939  } 

940  numberNodes2_++; 

941  nodeOnTree=true; // back on tree 

942  deleteNode = false ; 

943  # ifdef CHECK_NODE 

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

945  nodeInfo>numberBranchesLeft(), 

946  nodeInfo>numberPointingToThis()) ; 

947  # endif 

948  } 

949  else 

950  { deleteNode = true ; } 

951  

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

953  /* 

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

955  unless the current feasible solution region includes the optimal solution 

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

957  optimal path. 

958  */ 

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

960  if (debugger) 

961  { if(debugger>onOptimalPath(*solver_)) { 

962  onOptimalPath=true; 

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

964  } else { 

965  printf("Not on optimal path\n") ; } 

966  } 

967  } 

968  /* 

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

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

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

972  */ 

973  phase_=2; 

974  cuts = OsiCuts() ; 

975  currentNumberCuts = solver_>getNumRows()numberRowsAtContinuous_ ; 

976  int saveNumber = numberIterations_; 

977  feasible = solveWithCuts(cuts,maximumCutPasses_,node, 

978  numberOldActiveCuts,numberNewCuts, 

979  maximumWhich,whichGenerator) ; 

980  if (statistics) { 

981  assert (numberNodes2_); 

982  assert (statistics[numberNodes2_1]); 

983  assert (statistics[numberNodes2_1]>node()==numberNodes2_1); 

984  statistics[numberNodes2_1]>endOfBranch(numberIterations_saveNumber, 

985  feasible ? solver_>getObjValue() 

986  : COIN_DBL_MAX); 

987  } 

988  /* 

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

990  */ 

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

992  if (numberNodes_ < intParam_[CbcMaxNumNode] && 

993  numberSolutions_ < intParam_[CbcMaxNumSol] && 

994  totalTime < dblParam_[CbcMaximumSeconds] && 

995  !stoppedOnGap&&!eventHappened) 

996  { 

997  /* 

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

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

1000  monotone objects. 

1001  

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

1003  created back in solveWithCuts. addCuts() will also deal with the cut 

1004  reference counts. 

1005  

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

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

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

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

1010  */ 

1011  if (onOptimalPath) 

1012  assert (feasible); 

1013  if (feasible) 

1014  { newNode = new CbcNode ; 

1015  newNode>setObjectiveValue(direction*solver_>getObjValue()) ; 

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

1017  anyAction=2; 

1018  anyAction =1 ; 

1019  resolved = false ; 

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

1021  anyAction=2; 

1022  // only allow twenty passes 

1023  int numberPassesLeft=20; 

1024  while (anyAction == 1) 

1025  { 

1026  if (numberBeforeTrust_<=0 ) { 

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

1028  } else { 

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

1030  if (anyAction==3) 

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

1032  } 

1033  if (onOptimalPath) 

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

1035  numberPassesLeft; 

1036  if (numberPassesLeft<=1) { 

1037  if (!numberLongStrong) 

1038  messageHandler()>message(CBC_WARNING_STRONG, 

1039  messages()) << CoinMessageEol ; 

1040  numberLongStrong++; 

1041  } 

1042  if (anyAction == 1) 

1043  { 

1044  // can do quick optimality check 

1045  int easy=2; 

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

1047  feasible = resolve() ; 

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

1049  resolved = true ; 

1050  if (feasible) 

1051  { newNode>setObjectiveValue(direction* 

1052  solver_>getObjValue()) ; 

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

1054  anyAction=2; 

1055  } 

1056  else 

1057  { anyAction = 2 ; } } } 

1058  if (anyAction >= 0) 

1059  { if (resolved) 

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

1061  takeOffCuts(cuts,whichGenerator,numberOldActiveCuts, 

1062  numberNewCuts,needValidSolution) ; 

1063  # ifdef CHECK_CUT_COUNTS 

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

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

1066  numberRowsAtContinuous_+numberNewCuts+ 

1067  numberOldActiveCuts) ; 

1068  const CoinWarmStartBasis* debugws = 

1069  dynamic_cast<const CoinWarmStartBasis*> 

1070  (solver_>getWarmStart()) ; 

1071  debugws>print() ; 

1072  delete debugws ; } 

1073  # endif 

1074  } 

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

1076  numberOldActiveCuts,numberNewCuts) ; 

1077  if (newNode>numberUnsatisfied()) 

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

1079  whichGenerator) ; } } 

1080  else 

1081  { anyAction = 2 ; } 

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

1083  if ( anyAction >=0 ) { 

1084  assert (newNode); 

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

1086  anyAction = 2; // say bad after all 

1087  } 

1088  /* 

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

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

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

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

1093  */ 

1094  if (anyAction == 2) 

1095  { delete newNode ; 

1096  newNode = NULL ; 

1097  // switch off any hot start 

1098  hotstartStrategy_=0; 

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

1100  { if (addedCuts_[i]) 

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

1102  delete addedCuts_[i] ; } } } 

1103  else 

1104  { nodeInfo>increment() ; } 

1105  /* 

1106  At this point, there are three possibilities: 

1107  * We have a live node (variable() >= 0) which will require further 

1108  branching to resolve. Before we push it onto the search tree, try for 

1109  a heuristic solution. 

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

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

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

1113  deleted, and newNode is null. 

1114  

1115  TODO: (lh) Now I'm more confused. I thought that the call to addCuts() above 

1116  took care of incrementing the reference counts for cuts at newNode. 

1117  Clearly I need to look more carefully. 

1118  */ 

1119  if (eventHandler) { 

1120  if (!eventHandler>event(ClpEventHandler::node)) { 

1121  eventHappened=true; // exit 

1122  } 

1123  } 

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

1125  if (statistics) { 

1126  assert (numberNodes2_); 

1127  assert (statistics[numberNodes2_1]); 

1128  assert (statistics[numberNodes2_1]>node()==numberNodes2_1); 

1129  if (newNode) 

1130  statistics[numberNodes2_1]>updateInfeasibility(newNode>numberUnsatisfied()); 

1131  else 

1132  statistics[numberNodes2_1]>sayInfeasible(); 

1133  } 

1134  if (newNode) 

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

1136  { handler_>message(CBC_BRANCH,messages_) 

1137  << numberNodes_<< newNode>objectiveValue() 

1138  << newNode>numberUnsatisfied()<< newNode>depth() 

1139  << CoinMessageEol ; 

1140  // Increment cut counts (taking off current) 

1141  int numberLeft = newNode>numberBranches() ; 

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

1143  { if (addedCuts_[i]) 

1144  { 

1145  # ifdef CHECK_CUT_COUNTS 

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

1147  numberLeft1) ; 

1148  # endif 

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

1150  

1151  double estValue = newNode>guessedObjectiveValue() ; 

1152  int found = 1 ; 

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

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

1155  double * newSolution = new double [numberColumns] ; 

1156  double heurValue = getCutoff() ; 

1157  int iHeur ; 

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

1159  { double saveValue = heurValue ; 

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

1161  if (ifSol > 0) { 

1162  // new solution found 

1163  found = iHeur ; 

1164  incrementUsed(newSolution); 

1165  } 

1166  else 

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

1168  { estValue = CoinMin(heurValue,estValue) ; 

1169  heurValue = saveValue ; } } 

1170  if (found >= 0) 

1171  { setBestSolution(CBC_ROUNDING,heurValue,newSolution) ; } 

1172  delete [] newSolution ; 

1173  newNode>setGuessedObjectiveValue(estValue) ; 

1174  tree_>push(newNode) ; 

1175  if (statistics) { 

1176  if (numberNodes2_==maximumStatistics) { 

1177  maximumStatistics = 2*maximumStatistics; 

1178  CbcStatistics ** temp = new CbcStatistics * [maximumStatistics]; 

1179  memset(temp,0,maximumStatistics*sizeof(CbcStatistics *)); 

1180  memcpy(temp,statistics,numberNodes2_*sizeof(CbcStatistics *)); 

1181  delete [] statistics; 

1182  statistics=temp; 

1183  } 

1184  assert (!statistics[numberNodes2_]); 

1185  statistics[numberNodes2_]=new CbcStatistics(newNode); 

1186  } 

1187  numberNodes2_++; 

1188  # ifdef CHECK_NODE 

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

1190  # endif 

1191  } 

1192  else 

1193  { for (i = 0 ; i < currentNumberCuts_ ; i++) 

1194  { if (addedCuts_[i]) 

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

1196  delete addedCuts_[i] ; } } 

1197  double objectiveValue = newNode>objectiveValue(); 

1198  setBestSolution(CBC_SOLUTION,objectiveValue, 

1199  solver_>getColSolution()) ; 

1200  incrementUsed(solver_>getColSolution()); 

1201  assert(nodeInfo>numberPointingToThis() <= 2) ; 

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

1203  nodeInfo>increment(); 

1204  delete newNode ; 

1205  nodeInfo>decrement() ; } } 

1206  /* 

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

1208  set. 

1209  */ 

1210  if (deleteNode) 

1211  delete node ; } 

1212  /* 

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

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

1215  and break out of the node processing loop. 

1216  */ 

1217  else 

1218  { 

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

1220  delete nextRowCut_; 

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

1222  if (!nodeOnTree&&!stoppedOnGap&&node!=rootNode) 

1223  delete node; 

1224  if (stoppedOnGap) 

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

1226  << bestObjective_bestPossibleObjective_ 

1227  << dblParam_[CbcAllowableGap] 

1228  << dblParam_[CbcAllowableFractionGap]*100.0 

1229  << CoinMessageEol ; 

1230  secondaryStatus_ = 2; 

1231  status_ = 0 ; } 

1232  else 

1233  if (isNodeLimitReached()) 

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

1235  secondaryStatus_ = 3; 

1236  status_ = 1 ; } 

1237  else 

1238  if (totalTime >= dblParam_[CbcMaximumSeconds]) 

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

1240  secondaryStatus_ = 4; 

1241  status_ = 1 ; } 

1242  else 

1243  if (eventHappened) 

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

1245  secondaryStatus_ = 5; 

1246  status_ = 5 ; } 

1247  else 

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

1249  secondaryStatus_ = 6; 

1250  status_ = 1 ; } 

1251  break ; } 

1252  /* 

1253  Delete cuts to get back to the original system. 

1254  

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

1256  to this code block also performs this action. 

1257  */ 

1258  int numberToDelete = getNumRows()numberRowsAtContinuous_ ; 

1259  if (numberToDelete) 

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

1261  int i ; 

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

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

1264  solver_>deleteRows(numberToDelete,delRows) ; 

1265  delete [] delRows ; } } 

1266  /* 

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

1268  */ 

1269  else 

1270  { delete node ; } } 

1271  /* 

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

1273  we hit some limit on evaluation. 

1274  

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

1276  */ 

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

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

1279  tree_>endSearch(); 

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

1281  if ( numberStoppedSubTrees_) 

1282  status_=1; 

1283  if (!status_) { 

1284  bestPossibleObjective_=bestObjective_; 

1285  handler_>message(CBC_END_GOOD,messages_) 

1286  << bestObjective_ << numberIterations_ << numberNodes_ 

1287  << CoinMessageEol ; 

1288  } else { 

1289  handler_>message(CBC_END,messages_) 

1290  << bestObjective_ <<bestPossibleObjective_ 

1291  << numberIterations_ << numberNodes_ 

1292  << CoinMessageEol ; 

1293  } 

1294  if (statistics) { 

1295  // report in some way 

1296  int * lookup = new int[numberObjects_]; 

1297  int i; 

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

1299  lookup[i]=1; 

1300  bool goodIds=true; 

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

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

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

1304  if (iColumn<0) 

1305  iColumn = id+numberColumns; 

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

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

1308  lookup[id]=iColumn; 

1309  } else { 

1310  goodIds=false; 

1311  break; 

1312  } 

1313  } else { 

1314  goodIds=false; 

1315  break; 

1316  } 

1317  } 

1318  if (!goodIds) { 

1319  delete [] lookup; 

1320  lookup=NULL; 

1321  } 

1322  if (doStatistics==3) { 

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

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

1325  statistics[i]>print(lookup); 

1326  } 

1327  } 

1328  if (doStatistics>1) { 

1329  // Find last solution 

1330  int k; 

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

1332  if (statistics[k]>endingObjective()!=COIN_DBL_MAX&& 

1333  !statistics[k]>endingInfeasibility()) 

1334  break; 

1335  } 

1336  if (k>=0) { 

1337  int depth=statistics[k]>depth(); 

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

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

1340  which[i]=k; 

1341  k=statistics[k]>parentNode(); 

1342  } 

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

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

1345  statistics[which[i]]>print(lookup); 

1346  } 

1347  delete [] which; 

1348  } 

1349  } 

1350  // now summary 

1351  int maxDepth=0; 

1352  double averageSolutionDepth=0.0; 

1353  int numberSolutions=0; 

1354  double averageCutoffDepth=0.0; 

1355  double averageSolvedDepth=0.0; 

1356  int numberCutoff=0; 

1357  int numberDown=0; 

1358  int numberFirstDown=0; 

1359  double averageInfDown=0.0; 

1360  double averageObjDown=0.0; 

1361  int numberCutoffDown=0; 

1362  int numberUp=0; 

1363  int numberFirstUp=0; 

1364  double averageInfUp=0.0; 

1365  double averageObjUp=0.0; 

1366  int numberCutoffUp=0; 

1367  double averageNumberIterations1=0.0; 

1368  double averageValue=0.0; 

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

1370  int depth = statistics[i]>depth(); 

1371  int way = statistics[i]>way(); 

1372  double value = statistics[i]>value(); 

1373  double startingObjective = statistics[i]>startingObjective(); 

1374  int startingInfeasibility = statistics[i]>startingInfeasibility(); 

1375  double endingObjective = statistics[i]>endingObjective(); 

1376  int endingInfeasibility = statistics[i]>endingInfeasibility(); 

1377  maxDepth = CoinMax(depth,maxDepth); 

1378  // Only for completed 

1379  averageNumberIterations1 += statistics[i]>numberIterations(); 

1380  averageValue += value; 

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

1382  numberSolutions++; 

1383  averageSolutionDepth += depth; 

1384  } 

1385  if (endingObjective==COIN_DBL_MAX) { 

1386  numberCutoff++; 

1387  averageCutoffDepth += depth; 

1388  if (way<0) { 

1389  numberDown++; 

1390  numberCutoffDown++; 

1391  if (way==1) 

1392  numberFirstDown++; 

1393  } else { 

1394  numberUp++; 

1395  numberCutoffUp++; 

1396  if (way==1) 

1397  numberFirstUp++; 

1398  } 

1399  } else { 

1400  averageSolvedDepth += depth; 

1401  if (way<0) { 

1402  numberDown++; 

1403  averageInfDown += startingInfeasibilityendingInfeasibility; 

1404  averageObjDown += endingObjectivestartingObjective; 

1405  if (way==1) 

1406  numberFirstDown++; 

1407  } else { 

1408  numberUp++; 

1409  averageInfUp += startingInfeasibilityendingInfeasibility; 

1410  averageObjUp += endingObjectivestartingObjective; 

1411  if (way==1) 

1412  numberFirstUp++; 

1413  } 

1414  } 

1415  } 

1416  // Now print 

1417  if (numberSolutions) 

1418  averageSolutionDepth /= (double) numberSolutions; 

1419  int numberSolved = numberNodes2_numberCutoff; 

1420  double averageNumberIterations2=numberIterations_averageNumberIterations1; 

1421  if(numberCutoff) { 

1422  averageCutoffDepth /= (double) numberCutoff; 

1423  averageNumberIterations2 /= (double) numberCutoff; 

1424  } 

1425  if (numberNodes2_) 

1426  averageValue /= (double) numberNodes2_; 

1427  if (numberSolved) { 

1428  averageNumberIterations1 /= (double) numberSolved; 

1429  averageSolvedDepth /= (double) numberSolved; 

1430  } 

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

1432  numberSolutions,averageSolutionDepth); 

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

1434  averageValue); 

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

1436  numberCutoff,averageCutoffDepth,averageNumberIterations2); 

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

1438  numberSolved,averageSolvedDepth,averageNumberIterations1); 

1439  if (numberDown) { 

1440  averageInfDown /= (double) numberDown; 

1441  averageObjDown /= (double) numberDown; 

1442  } 

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

1444  numberDown,numberFirstDown,numberDownnumberFirstDown,numberCutoffDown, 

1445  averageInfDown,averageObjDown); 

1446  if (numberUp) { 

1447  averageInfUp /= (double) numberUp; 

1448  averageObjUp /= (double) numberUp; 

1449  } 

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

1451  numberUp,numberFirstUp,numberUpnumberFirstUp,numberCutoffUp, 

1452  averageInfUp,averageObjUp); 

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

1454  delete statistics[i]; 

1455  delete [] statistics; 

1456  delete [] lookup; 

1457  } 

1458  /* 

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

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

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

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

1463  

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

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

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

1467  tell us more. 

1468  

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

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

1471  public methods. 

1472  */ 

1473  if (bestSolution_) 

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

1475  phase_=5; 

1476  setBestSolution(CBC_SOLUTION,bestObjective_,bestSolution_,true) ; 

1477  continuousSolver_>resolve() ; 

1478  if (!continuousSolver_>isProvenOptimal()) 

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

1480  continuousSolver_>initialSolve() ; } 

1481  delete solver_ ; 

1482  solver_ = continuousSolver_ ; 

1483  continuousSolver_ = NULL ; } 

1484  /* 

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

1486  */ 

1487  delete lastws ; 

1488  delete [] whichGenerator ; 

1489  delete [] lowerBefore ; 

1490  delete [] upperBefore ; 

1491  delete [] walkback_ ; 

1492  walkback_ = NULL ; 

1493  delete [] addedCuts_ ; 

1494  addedCuts_ = NULL ; 

1495  if (continuousSolver_) 

1496  { delete continuousSolver_ ; 

1497  continuousSolver_ = NULL ; } 

1498  /* 

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

1500  */ 

1501  globalCuts_= OsiCuts() ; 

1502  return ; } 

1503  

1504  

1505  

1506  // Solve the initial LP relaxation 

1507  void 

1508  CbcModel::initialSolve() 

1509  { 

1510  assert (solver_); 

1511  solver_>initialSolve(); 

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

1513  status_=1; 

1514  secondaryStatus_ = 1; 

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

1516  delete [] continuousSolution_; 

1517  continuousSolution_ = CoinCopyOfArray(solver_>getColSolution(), 

1518  solver_>getNumCols()); 

1519  } 

1520  

1521  /*! \brief Get an empty basis object 

1522  

1523  Return an empty CoinWarmStartBasis object with the requested capacity, 

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

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

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

1527  result. 

1528  */ 

1529  

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

1531  

1532  { CoinWarmStartBasis *emptyBasis ; 

1533  /* 

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

1535  */ 

1536  if (emptyWarmStart_ == 0) 

1537  { if (solver_ == 0) 

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

1539  "getEmptyBasis","CbcModel") ; } 

1540  emptyBasis = 

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

1542  if (emptyBasis == 0) 

1543  { throw CoinError( 

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

1545  "getEmptyBasis","CbcModel") ; } 

1546  emptyBasis>setSize(0,0) ; 

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

1548  /* 

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

1550  */ 

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

1552  assert(emptyBasis) ; 

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

1554  

1555  return (emptyBasis) ; } 

1556  

1557  

1558  /** Default Constructor 

1559  

1560  Creates an empty model without an associated solver. 

1561  */ 

1562  CbcModel::CbcModel() 

1563  

1564  : 

1565  solver_(NULL), 

1566  ourSolver_(true), 

1567  continuousSolver_(NULL), 

1568  defaultHandler_(true), 

1569  emptyWarmStart_(NULL), 

1570  basis_(NULL), 

1571  bestObjective_(COIN_DBL_MAX), 

1572  bestPossibleObjective_(COIN_DBL_MAX), 

1573  sumChangeObjective1_(0.0), 

1574  sumChangeObjective2_(0.0), 

1575  bestSolution_(NULL), 

1576  currentSolution_(NULL), 

1577  testSolution_(NULL), 

1578  minimumDrop_(1.0e4), 

1579  numberSolutions_(0), 

1580  stateOfSearch_(0), 

1581  hotstartStrategy_(0), 

1582  numberHeuristicSolutions_(0), 

1583  numberNodes_(0), 

1584  numberNodes2_(0), 

1585  numberIterations_(0), 

1586  status_(1), 

1587  secondaryStatus_(1), 

1588  numberIntegers_(0), 

1589  numberRowsAtContinuous_(0), 

1590  maximumNumberCuts_(0), 

1591  phase_(0), 

1592  currentNumberCuts_(0), 

1593  maximumDepth_(0), 

1594  walkback_(NULL), 

1595  addedCuts_(NULL), 

1596  nextRowCut_(NULL), 

1597  currentNode_(NULL), 

1598  integerVariable_(NULL), 

1599  continuousSolution_(NULL), 

1600  usedInSolution_(NULL), 

1601  specialOptions_(0), 

1602  subTreeModel_(NULL), 

1603  numberStoppedSubTrees_(0), 

1604  presolve_(0), 

1605  numberStrong_(5), 

1606  numberBeforeTrust_(0), 

1607  numberPenalties_(20), 

1608  penaltyScaleFactor_(3.0), 

1609  numberInfeasibleNodes_(0), 

1610  problemType_(0), 

1611  printFrequency_(0), 

1612  numberCutGenerators_(0), 

1613  generator_(NULL), 

1614  virginGenerator_(NULL), 

1615  numberHeuristics_(0), 

1616  heuristic_(NULL), 

1617  numberObjects_(0), 

1618  object_(NULL), 

1619  originalColumns_(NULL), 

1620  howOftenGlobalScan_(1), 

1621  numberGlobalViolations_(0), 

1622  continuousObjective_(COIN_DBL_MAX), 

1623  originalContinuousObjective_(COIN_DBL_MAX), 

1624  continuousInfeasibilities_(INT_MAX), 

1625  maximumCutPassesAtRoot_(20), 

1626  maximumCutPasses_(10), 

1627  resolveAfterTakeOffCuts_(true) 

1628  { 

1629  intParam_[CbcMaxNumNode] = 99999999; 

1630  intParam_[CbcMaxNumSol] = 9999999; 

1631  intParam_[CbcFathomDiscipline] = 0; 

1632  

1633  dblParam_[CbcIntegerTolerance] = 1e6; 

1634  dblParam_[CbcInfeasibilityWeight] = 0.0; 

1635  dblParam_[CbcCutoffIncrement] = 1e5; 

1636  dblParam_[CbcAllowableGap] = 1.0e10; 

1637  dblParam_[CbcAllowableFractionGap] = 0.0; 

1638  dblParam_[CbcMaximumSeconds] = 1.0e100; 

1639  dblParam_[CbcStartSeconds] = 0.0; 

1640  nodeCompare_=new CbcCompareDefault();; 

1641  tree_= new CbcTree(); 

1642  branchingMethod_=NULL; 

1643  strategy_=NULL; 

1644  parentModel_=NULL; 

1645  appData_=NULL; 

1646  handler_ = new CoinMessageHandler(); 

1647  handler_>setLogLevel(2); 

1648  messages_ = CbcMessage(); 

1649  } 

1650  

1651  /** Constructor from solver. 

1652  

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

1654  */ 

1655  

1656  CbcModel::CbcModel(const OsiSolverInterface &rhs) 

1657  : 

1658  continuousSolver_(NULL), 

1659  defaultHandler_(true), 

1660  emptyWarmStart_(NULL), 

1661  basis_(NULL) , 

1662  bestObjective_(COIN_DBL_MAX), 

1663  bestPossibleObjective_(COIN_DBL_MAX), 

1664  sumChangeObjective1_(0.0), 

1665  sumChangeObjective2_(0.0), 

1666  minimumDrop_(1.0e4), 

1667  numberSolutions_(0), 

1668  stateOfSearch_(0), 

1669  hotstartStrategy_(0), 

1670  numberHeuristicSolutions_(0), 

1671  numberNodes_(0), 

1672  numberNodes2_(0), 

1673  numberIterations_(0), 

1674  status_(1), 

1675  secondaryStatus_(1), 

1676  numberRowsAtContinuous_(0), 

1677  maximumNumberCuts_(0), 

1678  phase_(0), 

1679  currentNumberCuts_(0), 

1680  maximumDepth_(0), 

1681  walkback_(NULL), 

1682  addedCuts_(NULL), 

1683  nextRowCut_(NULL), 

1684  currentNode_(NULL), 

1685  specialOptions_(0), 

1686  subTreeModel_(NULL), 

1687  numberStoppedSubTrees_(0), 

1688  presolve_(0), 

1689  numberStrong_(5), 

1690  numberBeforeTrust_(0), 

1691  numberPenalties_(20), 

1692  penaltyScaleFactor_(3.0), 

1693  numberInfeasibleNodes_(0), 

1694  problemType_(0), 

1695  printFrequency_(0), 

1696  numberCutGenerators_(0), 

1697  generator_(NULL), 

1698  virginGenerator_(NULL), 

1699  numberHeuristics_(0), 

1700  heuristic_(NULL), 

1701  numberObjects_(0), 

1702  object_(NULL), 

1703  originalColumns_(NULL), 

1704  howOftenGlobalScan_(1), 

1705  numberGlobalViolations_(0), 

1706  continuousObjective_(COIN_DBL_MAX), 

1707  originalContinuousObjective_(COIN_DBL_MAX), 

1708  continuousInfeasibilities_(INT_MAX), 

1709  maximumCutPassesAtRoot_(20), 

1710  maximumCutPasses_(10), 

1711  resolveAfterTakeOffCuts_(true) 

1712  { 

1713  intParam_[CbcMaxNumNode] = 99999999; 

1714  intParam_[CbcMaxNumSol] = 9999999; 

1715  intParam_[CbcFathomDiscipline] = 0; 

1716  

1717  dblParam_[CbcIntegerTolerance] = 1e6; 

1718  dblParam_[CbcInfeasibilityWeight] = 0.0; 

1719  dblParam_[CbcCutoffIncrement] = 1e5; 

1720  dblParam_[CbcAllowableGap] = 1.0e10; 

1721  dblParam_[CbcAllowableFractionGap] = 0.0; 

1722  dblParam_[CbcMaximumSeconds] = 1.0e100; 

1723  dblParam_[CbcStartSeconds] = 0.0; 

1724  

1725  nodeCompare_=new CbcCompareDefault();; 

1726  tree_= new CbcTree(); 

1727  branchingMethod_=NULL; 

1728  strategy_=NULL; 

1729  parentModel_=NULL; 

1730  appData_=NULL; 

1731  handler_ = new CoinMessageHandler(); 

1732  handler_>setLogLevel(2); 

1733  messages_ = CbcMessage(); 

1734  solver_ = rhs.clone(); 

1735  ourSolver_ = true ; 

1736  

1737  // Initialize solution and integer variable vectors 

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

1739  numberIntegers_=0; 

1740  int numberColumns = solver_>getNumCols(); 

1741  int iColumn; 

1742  if (numberColumns) { 

1743  // Space for current solution 

1744  currentSolution_ = new double[numberColumns]; 

1745  continuousSolution_ = new double[numberColumns]; 

1746  usedInSolution_ = new int[numberColumns]; 

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

1748  if( solver_>isInteger(iColumn)) 

1749  numberIntegers_++; 

1750  } 

1751  } else { 

1752  // empty model 

1753  currentSolution_=NULL; 

1754  continuousSolution_=NULL; 

1755  usedInSolution_=NULL; 

1756  } 

1757  testSolution_=currentSolution_; 

1758  if (numberIntegers_) { 

1759  integerVariable_ = new int [numberIntegers_]; 

1760  numberIntegers_=0; 

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

1762  if( solver_>isInteger(iColumn)) 

1763  integerVariable_[numberIntegers_++]=iColumn; 

1764  } 

1765  } else { 

1766  integerVariable_ = NULL; 

1767  } 

1768  } 

1769  

1770  /* 

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

1772  

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

1774  

1775  Assuming ownership matches usage in OsiSolverInterface 

1776  (cf. assignProblem, loadProblem). 

1777  

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

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

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

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

1782  here, maybe? 

1783  

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

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

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

1787  */ 

1788  void 

1789  CbcModel::assignSolver(OsiSolverInterface *&solver) 

1790  

1791  { 

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

1793  if (solver_) 

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

1795  

1796  if (ourSolver_) delete solver_ ; 

1797  solver_ = solver; 

1798  solver = NULL ; 

1799  ourSolver_ = true ; 

1800  /* 

1801  Basis information is solverspecific. 

1802  */ 

1803  if (basis_) 

1804  { delete basis_ ; 

1805  basis_ = 0 ; } 

1806  if (emptyWarmStart_) 

1807  { delete emptyWarmStart_ ; 

1808  emptyWarmStart_ = 0 ; } 

1809  /* 

1810  Initialize integer variable vector. 

1811  */ 

1812  numberIntegers_=0; 

1813  int numberColumns = solver_>getNumCols(); 

1814  int iColumn; 

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

1816  if( solver_>isInteger(iColumn)) 

1817  numberIntegers_++; 

1818  } 

1819  if (numberIntegers_) { 

1820  delete [] integerVariable_; 

1821  integerVariable_ = new int [numberIntegers_]; 

1822  numberIntegers_=0; 

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

1824  if( solver_>isInteger(iColumn)) 

1825  integerVariable_[numberIntegers_++]=iColumn; 

1826  } 

1827  } else { 

1828  integerVariable_ = NULL; 

1829  } 

1830  

1831  return ; 

1832  } 

1833  

1834  // Copy constructor. 

1835  

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

1837  : 

1838  continuousSolver_(NULL), 

1839  defaultHandler_(rhs.defaultHandler_), 

1840  emptyWarmStart_(NULL), 

1841  basis_(NULL), 

1842  bestObjective_(rhs.bestObjective_), 

1843  bestPossibleObjective_(rhs.bestPossibleObjective_), 

1844  sumChangeObjective1_(rhs.sumChangeObjective1_), 

1845  sumChangeObjective2_(rhs.sumChangeObjective2_), 

1846  minimumDrop_(rhs.minimumDrop_), 

1847  numberSolutions_(rhs.numberSolutions_), 

1848  stateOfSearch_(rhs.stateOfSearch_), 

1849  hotstartStrategy_(rhs.hotstartStrategy_), 

1850  numberHeuristicSolutions_(rhs.numberHeuristicSolutions_), 

1851  numberNodes_(rhs.numberNodes_), 

1852  numberNodes2_(rhs.numberNodes2_), 

1853  numberIterations_(rhs.numberIterations_), 

1854  status_(rhs.status_), 

1855  secondaryStatus_(rhs.secondaryStatus_), 

1856  specialOptions_(rhs.specialOptions_), 

1857  subTreeModel_(rhs.subTreeModel_), 

1858  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_), 

1859  presolve_(rhs.presolve_), 

1860  numberStrong_(rhs.numberStrong_), 

1861  numberBeforeTrust_(rhs.numberBeforeTrust_), 

1862  numberPenalties_(rhs.numberPenalties_), 

1863  penaltyScaleFactor_(penaltyScaleFactor_), 

1864  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_), 

1865  problemType_(rhs.problemType_), 

1866  printFrequency_(rhs.printFrequency_), 

1867  howOftenGlobalScan_(rhs.howOftenGlobalScan_), 

1868  numberGlobalViolations_(rhs.numberGlobalViolations_), 

1869  continuousObjective_(rhs.continuousObjective_), 

1870  originalContinuousObjective_(rhs.originalContinuousObjective_), 

1871  continuousInfeasibilities_(rhs.continuousInfeasibilities_), 

1872  maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_), 

1873  maximumCutPasses_( rhs.maximumCutPasses_), 

1874  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_) 

1875  { 

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

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

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

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

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

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

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

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

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

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

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

1887  if (rhs.basis_) basis_ = 

1888  dynamic_cast<CoinWarmStartBasis *>(rhs.basis_>clone()) ; 

1889  if (defaultHandler_) { 

1890  handler_ = new CoinMessageHandler(); 

1891  handler_>setLogLevel(2); 

1892  } else { 

1893  handler_ = rhs.handler_; 

1894  } 

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

1896  numberCutGenerators_ = rhs.numberCutGenerators_; 

1897  if (numberCutGenerators_) { 

1898  generator_ = new CbcCutGenerator * [numberCutGenerators_]; 

1899  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_]; 

1900  int i; 

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

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

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

1904  } 

1905  } else { 

1906  generator_=NULL; 

1907  virginGenerator_=NULL; 

1908  } 

1909  if (!noTree) 

1910  globalCuts_ = rhs.globalCuts_; 

1911  numberHeuristics_ = rhs.numberHeuristics_; 

1912  if (numberHeuristics_) { 

1913  heuristic_ = new CbcHeuristic * [numberHeuristics_]; 

1914  int i; 

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

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

1917  } 

1918  } else { 

1919  heuristic_=NULL; 

1920  } 

1921  numberObjects_=rhs.numberObjects_; 

1922  if (numberObjects_) { 

1923  object_ = new CbcObject * [numberObjects_]; 

1924  int i; 

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

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

1927  } else { 

1928  object_=NULL; 

1929  } 

1930  if (!noTree!rhs.continuousSolver_) 

1931  solver_ = rhs.solver_>clone(); 

1932  else 

1933  solver_ = rhs.continuousSolver_>clone(); 

1934  if (rhs.originalColumns_) { 

1935  int numberColumns = solver_>getNumCols(); 

1936  originalColumns_= new int [numberColumns]; 

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

1938  } else { 

1939  originalColumns_=NULL; 

1940  } 

1941  nodeCompare_=rhs.nodeCompare_>clone(); 

1942  tree_= rhs.tree_>clone(); 

1943  branchingMethod_=rhs.branchingMethod_; 

1944  if (rhs.strategy_) 

1945  strategy_=rhs.strategy_>clone(); 

1946  else 

1947  strategy_=NULL; 

1948  parentModel_=rhs.parentModel_; 

1949  appData_=rhs.appData_; 

1950  messages_ = rhs.messages_; 

1951  ourSolver_ = true ; 

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

1953  numberIntegers_=rhs.numberIntegers_; 

1954  if (numberIntegers_) { 

1955  integerVariable_ = new int [numberIntegers_]; 

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

1957  } else { 

1958  integerVariable_ = NULL; 

1959  } 

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

1961  int numberColumns = solver_>getNumCols(); 

1962  bestSolution_ = new double[numberColumns]; 

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

1964  } else { 

1965  bestSolution_=NULL; 

1966  } 

1967  if (!noTree) { 

1968  int numberColumns = solver_>getNumCols(); 

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

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

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

1972  } else { 

1973  currentSolution_=NULL; 

1974  continuousSolution_=NULL; 

1975  usedInSolution_=NULL; 

1976  } 

1977  testSolution_=currentSolution_; 

1978  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_; 

1979  maximumNumberCuts_=rhs.maximumNumberCuts_; 

1980  phase_ = rhs.phase_; 

1981  currentNumberCuts_=rhs.currentNumberCuts_; 

1982  maximumDepth_= rhs.maximumDepth_; 

1983  if (noTree) { 

1984  bestObjective_ = COIN_DBL_MAX; 

1985  numberSolutions_ =0; 

1986  stateOfSearch_= 0; 

1987  numberHeuristicSolutions_=0; 

1988  numberNodes_=0; 

1989  numberNodes2_=0; 

1990  numberIterations_=0; 

1991  status_=0; 

1992  subTreeModel_=NULL; 

1993  numberStoppedSubTrees_=0; 

1994  continuousObjective_=COIN_DBL_MAX; 

1995  originalContinuousObjective_=COIN_DBL_MAX; 

1996  continuousInfeasibilities_=INT_MAX; 

1997  maximumNumberCuts_=0; 

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

1999  bestPossibleObjective_ = COIN_DBL_MAX; 

2000  } 

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

2002  if (maximumNumberCuts_) { 

2003  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

2004  } else { 

2005  addedCuts_ = NULL; 

2006  } 

2007  nextRowCut_ = NULL; 

2008  currentNode_ = NULL; 

2009  if (maximumDepth_) 

2010  walkback_ = new CbcNodeInfo * [maximumDepth_]; 

2011  else 

2012  walkback_ = NULL; 

2013  synchronizeModel(); 

2014  } 

2015  

2016  // Assignment operator 

2017  CbcModel & 

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

2019  { 

2020  if (this!=&rhs) { 

2021  

2022  if (defaultHandler_) 

2023  { delete handler_; 

2024  handler_ = NULL; } 

2025  defaultHandler_ = rhs.defaultHandler_; 

2026  if (defaultHandler_) 

2027  { handler_ = new CoinMessageHandler(); 

2028  handler_>setLogLevel(2); } 

2029  else 

2030  { handler_ = rhs.handler_; } 

2031  messages_ = rhs.messages_; 

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

2033  

2034  delete solver_; 

2035  if (rhs.solver_) 

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

2037  else 

2038  { solver_ = 0 ; } 

2039  ourSolver_ = true ; 

2040  delete continuousSolver_ ; 

2041  if (rhs.continuousSolver_) 

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

2043  else 

2044  { continuousSolver_ = 0 ; } 

2045  

2046  delete emptyWarmStart_ ; 

2047  if (rhs.emptyWarmStart_) 

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

2049  else 

2050  { emptyWarmStart_ = 0 ; } 

2051  delete basis_ ; 

2052  if (rhs.basis_) 

2053  { basis_ = dynamic_cast<CoinWarmStartBasis *>(rhs.basis_>clone()) ; } 

2054  else 

2055  { basis_ = 0 ; } 

2056  

2057  bestObjective_ = rhs.bestObjective_; 

2058  bestPossibleObjective_=rhs.bestPossibleObjective_; 

2059  sumChangeObjective1_=rhs.sumChangeObjective1_; 

2060  sumChangeObjective2_=rhs.sumChangeObjective2_; 

2061  delete [] bestSolution_; 

2062  if (rhs.bestSolution_) { 

2063  int numberColumns = rhs.getNumCols(); 

2064  bestSolution_ = new double[numberColumns]; 

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

2066  } else { 

2067  bestSolution_=NULL; 

2068  } 

2069  int numberColumns = solver_>getNumCols(); 

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

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

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

2073  testSolution_=currentSolution_; 

2074  minimumDrop_ = rhs.minimumDrop_; 

2075  numberSolutions_=rhs.numberSolutions_; 

2076  stateOfSearch_= rhs.stateOfSearch_; 

2077  hotstartStrategy_=rhs.hotstartStrategy_; 

2078  numberHeuristicSolutions_=rhs.numberHeuristicSolutions_; 

2079  numberNodes_ = rhs.numberNodes_; 

2080  numberNodes2_ = rhs.numberNodes2_; 

2081  numberIterations_ = rhs.numberIterations_; 

2082  status_ = rhs.status_; 

2083  secondaryStatus_ = rhs.secondaryStatus_; 

2084  specialOptions_ = rhs.specialOptions_; 

2085  subTreeModel_ = rhs.subTreeModel_; 

2086  numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_; 

2087  presolve_ = rhs.presolve_; 

2088  numberStrong_ = rhs.numberStrong_; 

2089  numberBeforeTrust_ = rhs.numberBeforeTrust_; 

2090  numberPenalties_ = rhs.numberPenalties_; 

2091  penaltyScaleFactor_ = penaltyScaleFactor_; 

2092  numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_; 

2093  problemType_ = rhs.problemType_; 

2094  printFrequency_ = rhs.printFrequency_; 

2095  howOftenGlobalScan_=rhs.howOftenGlobalScan_; 

2096  numberGlobalViolations_=rhs.numberGlobalViolations_; 

2097  continuousObjective_=rhs.continuousObjective_; 

2098  originalContinuousObjective_ = rhs.originalContinuousObjective_; 

2099  continuousInfeasibilities_ = rhs.continuousInfeasibilities_; 

2100  maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_; 

2101  maximumCutPasses_ = rhs.maximumCutPasses_; 

2102  resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_; 

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

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

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

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

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

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

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

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

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

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

2113  globalCuts_ = rhs.globalCuts_; 

2114  int i; 

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

2116  delete generator_[i]; 

2117  delete virginGenerator_[i]; 

2118  } 

2119  delete [] generator_; 

2120  delete [] virginGenerator_; 

2121  delete [] heuristic_; 

2122  numberCutGenerators_ = rhs.numberCutGenerators_; 

2123  if (numberCutGenerators_) { 

2124  generator_ = new CbcCutGenerator * [numberCutGenerators_]; 

2125  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_]; 

2126  int i; 

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

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

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

2130  } 

2131  } else { 

2132  generator_=NULL; 

2133  virginGenerator_=NULL; 

2134  } 

2135  numberHeuristics_ = rhs.numberHeuristics_; 

2136  if (numberHeuristics_) { 

2137  heuristic_ = new CbcHeuristic * [numberHeuristics_]; 

2138  memcpy(heuristic_,rhs.heuristic_, 

2139  numberHeuristics_*sizeof(CbcHeuristic *)); 

2140  } else { 

2141  heuristic_=NULL; 

2142  } 

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

2144  delete object_[i]; 

2145  delete [] object_; 

2146  numberObjects_=rhs.numberObjects_; 

2147  if (numberObjects_) { 

2148  object_ = new CbcObject * [numberObjects_]; 

2149  int i; 

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

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

2152  } else { 

2153  object_=NULL; 

2154  } 

2155  delete [] originalColumns_; 

2156  if (rhs.originalColumns_) { 

2157  int numberColumns = rhs.getNumCols(); 

2158  originalColumns_= new int [numberColumns]; 

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

2160  } else { 

2161  originalColumns_=NULL; 

2162  } 

2163  nodeCompare_=rhs.nodeCompare_>clone(); 

2164  delete tree_; 

2165  tree_= rhs.tree_>clone(); 

2166  branchingMethod_=rhs.branchingMethod_; 

2167  delete strategy_; 

2168  if (rhs.strategy_) 

2169  strategy_=rhs.strategy_>clone(); 

2170  else 

2171  strategy_=NULL; 

2172  parentModel_=rhs.parentModel_; 

2173  appData_=rhs.appData_; 

2174  

2175  delete [] integerVariable_; 

2176  numberIntegers_=rhs.numberIntegers_; 

2177  if (numberIntegers_) { 

2178  integerVariable_ = new int [numberIntegers_]; 

2179  memcpy(integerVariable_,rhs.integerVariable_, 

2180  numberIntegers_*sizeof(int)); 

2181  } else { 

2182  integerVariable_ = NULL; 

2183  } 

2184  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_; 

2185  maximumNumberCuts_=rhs.maximumNumberCuts_; 

2186  phase_ = rhs.phase_; 

2187  currentNumberCuts_=rhs.currentNumberCuts_; 

2188  maximumDepth_= rhs.maximumDepth_; 

2189  delete [] addedCuts_; 

2190  delete [] walkback_; 

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

2192  if (maximumNumberCuts_) { 

2193  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

2194  } else { 

2195  addedCuts_ = NULL; 

2196  } 

2197  nextRowCut_ = NULL; 

2198  currentNode_ = NULL; 

2199  if (maximumDepth_) 

2200  walkback_ = new CbcNodeInfo * [maximumDepth_]; 

2201  else 

2202  walkback_ = NULL; 

2203  synchronizeModel(); 

2204  } 

2205  return *this; 

2206  } 

2207  

2208  // Destructor 

2209  CbcModel::~CbcModel () 

2210  { 

2211  if (defaultHandler_) { 

2212  delete handler_; 

2213  handler_ = NULL; 

2214  } 

2215  delete tree_; 

2216  if (ourSolver_) delete solver_; 

2217  gutsOfDestructor(); 

2218  } 

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

2220  void 

2221  CbcModel::gutsOfDestructor() 

2222  { 

2223  delete emptyWarmStart_ ; 

2224  emptyWarmStart_ =NULL; 

2225  delete basis_ ; 

2226  basis_ =NULL; 

2227  delete continuousSolver_; 

2228  continuousSolver_=NULL; 

2229  delete [] bestSolution_; 

2230  bestSolution_=NULL; 

2231  delete [] currentSolution_; 

2232  currentSolution_=NULL; 

2233  delete [] continuousSolution_; 

2234  continuousSolution_=NULL; 

2235  delete [] usedInSolution_; 

2236  usedInSolution_ = NULL; 

2237  testSolution_=NULL; 

2238  delete [] integerVariable_; 

2239  integerVariable_=NULL; 

2240  int i; 

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

2242  delete generator_[i]; 

2243  delete virginGenerator_[i]; 

2244  } 

2245  delete [] generator_; 

2246  delete [] virginGenerator_; 

2247  generator_=NULL; 

2248  virginGenerator_=NULL; 

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

2250  delete heuristic_[i]; 

2251  delete [] heuristic_; 

2252  heuristic_=NULL; 

2253  delete nodeCompare_; 

2254  nodeCompare_=NULL; 

2255  delete [] addedCuts_; 

2256  addedCuts_=NULL; 

2257  nextRowCut_ = NULL; 

2258  currentNode_ = NULL; 

2259  delete [] walkback_; 

2260  walkback_=NULL; 

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

2262  delete object_[i]; 

2263  delete [] object_; 

2264  object_=NULL; 

2265  delete [] originalColumns_; 

2266  originalColumns_=NULL; 

2267  delete strategy_; 

2268  } 

2269  // Are there a numerical difficulties? 

2270  bool 

2271  CbcModel::isAbandoned() const 

2272  { 

2273  return status_ == 2; 

2274  } 

2275  // Is optimality proven? 

2276  bool 

2277  CbcModel::isProvenOptimal() const 

2278  { 

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

2280  return true; 

2281  else 

2282  return false; 

2283  } 

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

2285  bool 

2286  CbcModel::isProvenInfeasible() const 

2287  { 

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

2289  return true; 

2290  else 

2291  return false; 

2292  } 

2293  // Node limit reached? 

2294  bool 

2295  CbcModel::isNodeLimitReached() const 

2296  { 

2297  return numberNodes_ >= intParam_[CbcMaxNumNode]; 

2298  } 

2299  // Time limit reached? 

2300  bool 

2301  CbcModel::isSecondsLimitReached() const 

2302  { 

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

2304  return true; 

2305  else 

2306  return false; 

2307  } 

2308  // Solution limit reached? 

2309  bool 

2310  CbcModel::isSolutionLimitReached() const 

2311  { 

2312  return numberSolutions_ >= intParam_[CbcMaxNumSol]; 

2313  } 

2314  // Set language 

2315  void 

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

2317  { 

2318  messages_ = CbcMessage(language); 

2319  } 

2320  void 

2321  CbcModel::setNumberStrong(int number) 

2322  { 

2323  if (number<0) 

2324  numberStrong_=0; 

2325  else 

2326  numberStrong_=number; 

2327  } 

2328  void 

2329  CbcModel::setNumberBeforeTrust(int number) 

2330  { 

2331  if (number<=0) { 

2332  numberBeforeTrust_=0; 

2333  } else { 

2334  numberBeforeTrust_=number; 

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

2336  } 

2337  } 

2338  void 

2339  CbcModel::setNumberPenalties(int number) 

2340  { 

2341  if (number<=0) { 

2342  numberPenalties_=0; 

2343  } else { 

2344  numberPenalties_=number; 

2345  } 

2346  } 

2347  void 

2348  CbcModel::setPenaltyScaleFactor(double value) 

2349  { 

2350  if (value<=0) { 

2351  penaltyScaleFactor_=3.0; 

2352  } else { 

2353  penaltyScaleFactor_=value; 

2354  } 

2355  } 

2356  void 

2357  CbcModel::setHowOftenGlobalScan(int number) 

2358  { 

2359  if (number<1) 

2360  howOftenGlobalScan_=0; 

2361  else 

2362  howOftenGlobalScan_=number; 

2363  } 

2364  

2365  // Add one generator 

2366  void 

2367  CbcModel::addCutGenerator(CglCutGenerator * generator, 

2368  int howOften, const char * name, 

2369  bool normal, bool atSolution, 

2370  bool whenInfeasible,int howOftenInSub, 

2371  int whatDepth, int whatDepthInSub) 

2372  { 

2373  CbcCutGenerator ** temp = generator_; 

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

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

2376  delete[] temp ; 

2377  generator_[numberCutGenerators_]= 

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

2379  normal,atSolution,whenInfeasible,howOftenInSub, 

2380  whatDepth, whatDepthInSub); 

2381  // and before any cahnges 

2382  temp = virginGenerator_; 

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

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

2385  delete[] temp ; 

2386  virginGenerator_[numberCutGenerators_++]= 

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

2388  normal,atSolution,whenInfeasible,howOftenInSub, 

2389  whatDepth, whatDepthInSub); 

2390  

2391  } 

2392  // Add one heuristic 

2393  void 

2394  CbcModel::addHeuristic(CbcHeuristic * generator) 

2395  { 

2396  CbcHeuristic ** temp = heuristic_; 

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

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

2399  delete [] temp; 

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

2401  } 

2402  

2403  /* 

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

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

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

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

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

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

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

2411  applyToModel does the heavy lifting. 

2412  

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

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

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

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

2417  collection and model adjustment. 

2418  

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

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

2421  */ 

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

2423  { int i; 

2424  int nNode=0; 

2425  int numberColumns = getNumCols(); 

2426  CbcNodeInfo * nodeInfo = node>nodeInfo(); 

2427  

2428  /* 

2429  Remove all cuts from the constraint system. 

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

2431  the reference isn't clear to me). 

2432  */ 

2433  int currentNumberCuts = solver_>getNumRows()numberRowsAtContinuous_; 

2434  int *which = new int[currentNumberCuts]; 

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

2436  which[i] = i+numberRowsAtContinuous_; 

2437  solver_>deleteRows(currentNumberCuts,which); 

2438  delete [] which; 

2439  /* 

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

2441  cut count in currentNumberCuts. 

2442  

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

2444  old node (for cuts?) 

2445  */ 

2446  currentNumberCuts = 0; 

2447  while (nodeInfo) { 

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

2449  walkback_[nNode++]=nodeInfo; 

2450  currentNumberCuts += nodeInfo>numberCuts() ; 

2451  nodeInfo = nodeInfo>parent() ; 

2452  if (nNode==maximumDepth_) { 

2453  maximumDepth_ *= 2; 

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

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

2456  temp[i] = walkback_[i]; 

2457  delete [] walkback_; 

2458  walkback_ = temp; 

2459  } 

2460  } 

2461  /* 

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

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

2464  record those cuts in addedCuts_. 

2465  

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

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

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

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

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

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

2472  */ 

2473  currentNumberCuts_=currentNumberCuts; 

2474  if (currentNumberCuts >= maximumNumberCuts_) { 

2475  maximumNumberCuts_ = currentNumberCuts; 

2476  delete [] addedCuts_; 

2477  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

2478  } 

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

2480  /* 

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

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

2483  * lastws will be an appropriate basis for node; 

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

2485  node; and 

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

2487  constraint system at node. 

2488  applyToModel does all the heavy lifting. 

2489  */ 

2490  currentNumberCuts=0; 

2491  while (nNode) { 

2492  nNode; 

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

2494  } 

2495  } 

2496  

2497  /* 

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

2499  the cuts we've collected, add the ones that are tight and omit the ones that 

2500  are loose. If the node is infeasible, we just adjust the reference counts to 

2501  reflect that we're about to prune this node and its descendants. 

2502  

2503  The reason we need to pass in lastws is that OsiClp automagically corrects 

2504  the basis when it deletes constraints. So when all cuts are stripped within 

2505  addCuts1, we lose their basis entries, hence the ability to determine if 

2506  they are loose or tight. The question is whether we really need to pass in 

2507  a basis or if we can capture it here. I'm thinking we can capture it here 

2508  and pass it back out if required. 

2509  */ 

2510  int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws) 

2511  { 

2512  /* 

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

2514  comments there. 

2515  */ 

2516  addCuts1(node,lastws); 

2517  int i; 

2518  int numberColumns = getNumCols(); 

2519  CbcNodeInfo * nodeInfo = node>nodeInfo(); 

2520  double cutoff = getCutoff() ; 

2521  int currentNumberCuts=currentNumberCuts_; 

2522  /* 

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

2524  constraint system. 

2525  */ 

2526  if (node>objectiveValue() < cutoff) 

2527  { int numberToAdd = 0; 

2528  const OsiRowCut * * addCuts; 

2529  if (currentNumberCuts == 0) 

2530  addCuts = NULL; 

2531  else 

2532  addCuts = new const OsiRowCut * [currentNumberCuts]; 

2533  # ifdef CHECK_CUT_COUNTS 

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

2535  numberRowsAtContinuous_,currentNumberCuts); 

2536  lastws>print(); 

2537  # endif 

2538  /* 

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

2540  There are three steps: 

2541  1) Scan the basis. If the logical associated with the cut is basic, it's 

2542  loose and we drop it. The status of the logical for tight cuts is 

2543  written back into the status array, compressing as we go. 

2544  2) Resize the basis to fit the number of active cuts, stash a clone, and 

2545  install with a call to setWarmStart(). 

2546  3) Install the tight cuts into the constraint system (applyRowCuts). 

2547  

2548  TODO: After working through the code in createInfo, I'm more comfortable if 

2549  inactive cuts are retained in lastws. So, instead of cloning 

2550  lastws into basis_ after the compression loop, do it ahead of time 

2551  and then recover lastws from basis_ after the setWarmStart(). 

2552  (Minimal code change :). See CbcNode::createInfo for more. 

2553  */ 

2554  if (basis_) delete basis_ ; 

2555  basis_= dynamic_cast<CoinWarmStartBasis *>(lastws>clone()) ; 

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

2557  CoinWarmStartBasis::Status status = 

2558  lastws>getArtifStatus(i+numberRowsAtContinuous_); 

2559  if (status != CoinWarmStartBasis::basic&&addedCuts_[i]) { 

2560  # ifdef CHECK_CUT_COUNTS 

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

2562  numberRowsAtContinuous_+numberToAdd); 

2563  # endif 

2564  lastws>setArtifStatus(numberToAdd+numberRowsAtContinuous_,status); 

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

2566  } else { 

2567  # ifdef CHECK_CUT_COUNTS 

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

2569  # endif 

2570  addedCuts_[i]=NULL; 

2571  } 

2572  } 

2573  int numberRowsNow=numberRowsAtContinuous_+numberToAdd; 

2574  lastws>resize(numberRowsNow,numberColumns); 

2575  #ifdef FULL_DEBUG 

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

2577  numberRowsAtContinuous_,numberToAdd); 

2578  lastws>print(); 

2579  #endif 

2580  /* 

2581  Apply the cuts and set the basis in the solver. 

2582  */ 

2583  solver_>applyRowCuts(numberToAdd,addCuts); 

2584  solver_>setWarmStart(lastws); 

2585  /* 

2586  TODO: Undo the debugging change. Delete lastws and assign basis_. 

2587  */ 

2588  delete lastws ; 

2589  lastws = basis_ ; 

2590  basis_ = 0 ; 

2591  

2592  #if 0 

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

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

2595  node>objectiveValue(), 

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

2597  } 

2598  #endif 

2599  /* 

2600  Clean up and we're out of here. 

2601  */ 

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

2603  delete addCuts[i]; 

2604  delete [] addCuts; 

2605  numberNodes_++; 

2606  return 0; 

2607  } 

2608  /* 

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

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

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

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

2613  */ 

2614  else 

2615  { int i; 

2616  int numberLeft = nodeInfo>numberBranchesLeft(); 

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

2618  { if (addedCuts_[i]) 

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

2620  { delete addedCuts_[i]; 

2621  addedCuts_[i] = NULL; } } } 

2622  return 1 ; } 

2623  } 

2624  

2625  

2626  /* 

2627  Perform reduced cost fixing on integer variables. 

2628  

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

2630  down the current situation. 

2631  */ 

2632  

2633  void CbcModel::reducedCostFix () 

2634  

2635  { double cutoff = getCutoff() ; 

2636  double direction = solver_>getObjSense() ; 

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

2638  double integerTolerance = getDblParam(CbcIntegerTolerance) ; 

2639  

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

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

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

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

2644  

2645  int numberFixed = 0 ; 

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

2647  { int iColumn = integerVariable_[i] ; 

2648  double djValue = direction*reducedCost[iColumn] ; 

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

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

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

2652  numberFixed++ ; } 

2653  else 

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

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

2656  numberFixed++ ; } } } 

2657  

2658  return ; } 

2659  

2660  

2661  

2662  

2663  /** Solve the model using cuts 

2664  

2665  This version takes off redundant cuts from node. 

2666  Returns true if feasible. 

2667  

2668  \todo 

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

2670  relaxation and calling solveWithCuts? 

2671  

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

2673  */ 

2674  

2675  bool 

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

2677  int &numberOldActiveCuts, int &numberNewCuts, 

2678  int &maximumWhich, int *&whichGenerator) 

2679  /* 

2680  Parameters: 

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

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

2683  

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

2685  whichGenerator: (i/o) whichGenerator[i] is loaded with the index of the 

2686  generator that produced cuts[i]; reallocated as 

2687  required 

2688  numberOldActiveCuts: (o) the number of active cuts at this node from 

2689  previous rounds of cut generation 

2690  numberNewCuts: (o) the number of cuts produced in this round of cut 

2691  generation 

2692  maximumWhich: (i/o) capacity of whichGenerator; may be updated if 

2693  whichGenerator grows. 

2694  

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

2696  */ 

2697  

2698  

2699  { bool feasible = true ; 

2700  int lastNumberCuts = 0 ; 

2701  double lastObjective = 1.0e100 ; 

2702  int violated = 0 ; 

2703  int numberRowsAtStart = solver_>getNumRows() ; 

2704  int numberColumns = solver_>getNumCols() ; 

2705  

2706  numberOldActiveCuts = numberRowsAtStartnumberRowsAtContinuous_ ; 

2707  numberNewCuts = 0 ; 

2708  

2709  bool onOptimalPath = false ; 

2710  const OsiRowCutDebugger *debugger = NULL; 

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

2712  /* 

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

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

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

2716  */ 

2717  debugger = solver_>getRowCutDebugger() ; 

2718  if (debugger) 

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

2720  } 

2721  /* 

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

2723  after the debug stuff. 

2724  */ 

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

2726  if (node) 

2727  objectiveValue= node>objectiveValue(); 

2728  feasible = resolve() ; 

2729  

2730  // Update branching information if wanted 

2731  if(node &&branchingMethod_) 

2732  branchingMethod_>updateInformation(solver_,node); 

2733  

2734  #ifdef CBC_DEBUG 

2735  if (feasible) 

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

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

2738  solver_>getNumRows()) ; } 

2739  

2740  else 

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

2742  #endif 

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

2744  /* 

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

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

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

2748  */ 

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

2750  if (!feasible) { 

2751  solver_>writeMps("infeas"); 

2752  CoinWarmStartBasis *slack = 

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

2754  solver_>setWarmStart(slack); 

2755  delete slack ; 

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

2757  solver_>initialSolve(); 

2758  } 

2759  assert(feasible) ; 

2760  } 

2761  } 

2762  

2763  if (!feasible) { 

2764  numberInfeasibleNodes_++; 

2765  return (false) ; 

2766  } 

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

2768   objectiveValue ; 

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

2770  numberTries=0; // exit 

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

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

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

2774  // Return at once if numberTries zero 

2775  if (!numberTries) { 

2776  cuts=OsiCuts(); 

2777  numberNewCuts=0; 

2778  return true; 

2779  } 

2780  /* 

2781  Do reduced cost fixing, and then grab the primal solution and bounds vectors. 

2782  */ 

2783  reducedCostFix() ; 

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

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

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

2787  /* 

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

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

2790  the specified number of rounds. 

2791  */ 

2792  double minimumDrop = minimumDrop_ ; 

2793  if (numberTries<0) 

2794  { numberTries = numberTries ; 

2795  minimumDrop = 1.0 ; } 

2796  /* 

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

2798  up to do it. 

2799  */ 

2800  # define SCANCUTS 100 

2801  int *countColumnCuts = NULL ; 

2802  // Always accumulate row cut counts 

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

2804  memset(countRowCuts,0, 

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

2806  bool fullScan = false ; 

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

2808  { fullScan = true ; 

2809  countColumnCuts = new int[numberCutGenerators_+numberHeuristics_] ; 

2810  memset(countColumnCuts,0, 

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

2812  

2813  double direction = solver_>getObjSense() ; 

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

2815  

2816  currentPassNumber_ = 0 ; 

2817  double primalTolerance = 1.0e7 ; 

2818  /* 

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

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

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

2822  consider scanning globalCuts_ and activating any applicable cuts. 

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

2824  generator_ and heuristic_ arrays. Newly generated global cuts are 

2825  copied to globalCuts_ at this time. 

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

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

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

2829  does the necessary bookkeeping in the model. 

2830  */ 

2831  do 

2832  { currentPassNumber_++ ; 

2833  numberTries ; 

2834  OsiCuts theseCuts ; 

2835  /* 

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

2837  useful. 

2838  I can't see why this code 

2839  needs its own copy of the primal solution. Removed the dec'l. 

2840  */ 

2841  int numberViolated=0; 

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

2843  (numberNodes_%howOftenGlobalScan_) == 0) 

2844  { int numberCuts = globalCuts_.sizeColCuts() ; 

2845  int i; 

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

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

2848  if (thisCut>violated(solution)>primalTolerance) { 

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

2850  thisCut>violated(solution)) ; 

2851  whichGenerator[numberViolated++]=1; 

2852  theseCuts.insert(*thisCut) ; 

2853  } 

2854  } 

2855  numberCuts = globalCuts_.sizeRowCuts() ; 

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

2857  const OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ; 

2858  if (thisCut>violated(solution)>primalTolerance) { 

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

2860  // thisCut>violated(solution)) ; 

2861  whichGenerator[numberViolated++]=1; 

2862  theseCuts.insert(*thisCut) ; 

2863  } 

2864  } 

2865  numberGlobalViolations_+=numberViolated; 

2866  } 

2867  /* 

2868  Generate new cuts (global and/or local) and/or apply heuristics. If 

2869  CglProbing is used, then it should be first as it can fix continuous 

2870  variables. 

2871  

2872  At present, CglProbing is the only case where generateCuts will return 

2873  true. generateCuts actually modifies variable bounds in the solver when 

2874  CglProbing indicates that it can fix a variable. Reoptimisation is required 

2875  to take full advantage. 

2876  */ 

2877  if (nextRowCut_) { 

2878  // branch was a cut  add it 

2879  theseCuts.insert(*nextRowCut_); 

2880  //nextRowCut_>print(); 

2881  const OsiRowCut * cut=nextRowCut_; 

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

2883  double lb = cut>lb(); 

2884  double ub = cut>ub(); 

2885  int n=cut>row().getNumElements(); 

2886  const int * column = cut>row().getIndices(); 

2887  const double * element = cut>row().getElements(); 

2888  double sum=0.0; 

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

2890  int iColumn = column[i]; 

2891  double value = element[i]; 

2892  //if (solution[iColumn]>1.0e7) 

2893  //printf("value of %d is %g\n",iColumn,solution[iColumn]); 

2894  sum += value * solution[iColumn]; 

2895  } 

2896  delete nextRowCut_; 

2897  nextRowCut_=NULL; 

2898  if (handler_>logLevel()>1) 

2899  printf("applying branch cut, sum is %g, bounds %g %g\n",sum,lb,ub); 

2900  // set whichgenerator (also serves as marker to say don't delete0 

2901  whichGenerator[numberViolated++]=2; 

2902  } 

2903  double * newSolution = new double [numberColumns] ; 

2904  double heuristicValue = getCutoff() ; 

2905  int found = 1; // no solution found 

2906  

2907  for (int i = 0;i<numberCutGenerators_+numberHeuristics_;i++) { 

2908  int numberRowCutsBefore = theseCuts.sizeRowCuts() ; 

2909  int numberColumnCutsBefore = theseCuts.sizeColCuts() ; 

2910  if (i<numberCutGenerators_) { 

2911  if (generator_[i]>normal()) { 

2912  bool mustResolve = 

2913  generator_[i]>generateCuts(theseCuts,fullScan,node) ; 

2914  #ifdef CBC_DEBUG 

2915  { 

2916  int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 

2917  int k ; 

2918  for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 

2919  OsiRowCut thisCut = theseCuts.rowCut(k) ; 

2920  /* check size of elements. 

2921  We can allow smaller but this helps debug generators as it 

2922  is unsafe to have small elements */ 

2923  int n=thisCut.row().getNumElements(); 

2924  const int * column = thisCut.row().getIndices(); 

2925  const double * element = thisCut.row().getElements(); 

2926  //assert (n); 

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

2928  int iColumn = column[i]; 

2929  double value = element[i]; 

2930  assert(fabs(value)>1.0e12&&fabs(value)<1.0e20); 

2931  } 

2932  } 

2933  } 

2934  #endif 

2935  if (mustResolve) { 

2936  feasible = resolve() ; 

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

2938  debugger = solver_>getRowCutDebugger() ; 

2939  if (debugger) 

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

2941  else 

2942  onOptimalPath=false; 

2943  if (onOptimalPath && !solver_>isDualObjectiveLimitReached()) 

2944  assert(feasible) ; 

2945  } 

2946  if (!feasible) 

2947  break ; 

2948  } 

2949  } 

2950  } else { 

2951  // see if heuristic will do anything 

2952  double saveValue = heuristicValue ; 

2953  int ifSol = 

2954  heuristic_[inumberCutGenerators_]>solution(heuristicValue, 

2955  newSolution, 

2956  theseCuts) ; 

2957  if (ifSol>0) { 

2958  // better solution found 

2959  found = i ; 

2960  incrementUsed(newSolution); 

2961  } else if (ifSol<0) { 

2962  heuristicValue = saveValue ; 

2963  } 

2964  } 

2965  int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 

2966  int numberColumnCutsAfter = theseCuts.sizeColCuts() ; 

2967  

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

2969  if (onOptimalPath) { 

2970  int k ; 

2971  for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 

2972  OsiRowCut thisCut = theseCuts.rowCut(k) ; 

2973  if(debugger>invalidCut(thisCut)) { 

2974  solver_>writeMps("badCut"); 

2975  } 

2976  assert(!debugger>invalidCut(thisCut)) ; 

2977  } 

2978  } 

2979  } 

2980  

2981  /* 

2982  The cut generator/heuristic has done its thing, and maybe it generated some 

2983  cuts and/or a new solution. Do a bit of bookkeeping: load 

2984  whichGenerator[i] with the index of the generator responsible for a cut, 

2985  and place cuts flagged as global in the global cut pool for the model. 

2986  

2987  lastNumberCuts is the sum of cuts added in previous iterations; it's the 

2988  offset to the proper starting position in whichGenerator. 

2989  

2990  TODO: Why is whichGenerator increased to 2*maximumWhich when it grows? 

2991  */ 

2992  int numberBefore = 

2993  numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ; 

2994  int numberAfter = 

2995  numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ; 

2996  if (numberAfter > maximumWhich) { 

2997  maximumWhich = CoinMax(maximumWhich*2+100,numberAfter) ; 

2998  int * temp = new int[2*maximumWhich] ; 

2999  memcpy(temp,whichGenerator,numberBefore*sizeof(int)) ; 

3000  delete [] whichGenerator ; 

3001  whichGenerator = temp ; 

3002  } 

3003  int j ; 

3004  if (fullScan) { 

3005  // counts 

3006  countColumnCuts[i] += numberColumnCutsAfternumberColumnCutsBefore ; 

3007  } 

3008  countRowCuts[i] += numberRowCutsAfternumberRowCutsBefore ; 

3009  

3010  for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) { 

3011  whichGenerator[numberBefore++] = i ; 

3012  const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ; 

3013  if (thisCut>lb()>thisCut>ub()) 

3014  violated=2; // subproblem is infeasible 

3015  if (thisCut>globallyValid()) { 

3016  // add to global list 

3017  globalCuts_.insert(*thisCut) ; 

3018  } 

3019  } 

3020  for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) { 

3021  whichGenerator[numberBefore++] = i ; 

3022  const OsiColCut * thisCut = theseCuts.colCutPtr(j) ; 

3023  if (thisCut>globallyValid()) { 

3024  // add to global list 

3025  globalCuts_.insert(*thisCut) ; 

3026  } 

3027  } 

3028  } 

3029  // If at root node and first pass do heuristics without cuts 

3030  if (!numberNodes_&¤tPassNumber_==1) { 

3031  // Save number solutions 

3032  int saveNumberSolutions = numberSolutions_; 

3033  int saveNumberHeuristicSolutions = numberHeuristicSolutions_; 

3034  for (int i = 0;i<numberHeuristics_;i++) { 

3035  // see if heuristic will do anything 

3036  double saveValue = heuristicValue ; 

3037  int ifSol = heuristic_[i]>solution(heuristicValue, 

3038  newSolution); 

3039  if (ifSol>0) { 

3040  // better solution found 

3041  found = i ; 

3042  incrementUsed(newSolution); 

3043  // increment number of solutions so other heuristics can test 

3044  numberSolutions_++; 

3045  numberHeuristicSolutions_++; 

3046  } else { 

3047  heuristicValue = saveValue ; 

3048  } 

3049  } 

3050  // Restore number solutions 

3051  numberSolutions_ = saveNumberSolutions; 

3052  numberHeuristicSolutions_ = saveNumberHeuristicSolutions; 

3053  } 

3054  /* 

3055  End of the loop to exercise each generator/heuristic. 

3056  

3057  Did any of the heuristics turn up a new solution? Record it before we free 

3058  the vector. 

3059  */ 

3060  if (found >= 0) 

3061  { 

3062  phase_=4; 

3063  incrementUsed(newSolution); 

3064  setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; } 

3065  delete [] newSolution ; 

3066  

3067  #if 0 

3068  // switch on to get all cuts printed 

3069  theseCuts.printCuts() ; 

3070  #endif 

3071  int numberColumnCuts = theseCuts.sizeColCuts() ; 

3072  int numberRowCuts = theseCuts.sizeRowCuts() ; 

3073  if (violated>=0) 

3074  violated = numberRowCuts + numberColumnCuts ; 

3075  /* 

3076  Apply column cuts (aka bound tightening). This may be partially redundant 

3077  for column cuts returned by CglProbing, as generateCuts installs bounds 

3078  from CglProbing when it determines it can fix a variable. 

3079  

3080  TODO: Looks like the use of violated has evolved. The value set above is 

3081  completely ignored. All that's left is violated == 1 indicates some 

3082  cut is violated, violated == 2 indicates infeasibility. Only 

3083  infeasibility warrants exceptional action. 

3084  

3085  TODO: Strikes me that this code will fail to detect infeasibility, because 

3086  the breaks escape the inner loops but the outer loop keeps going. 

3087  Infeasibility in an early cut will be overwritten if a later cut is 

3088  merely violated. 

3089  */ 

3090  if (numberColumnCuts) { 

3091  

3092  #ifdef CBC_DEBUG 

3093  double * oldLower = new double [numberColumns] ; 

3094  double * oldUpper = new double [numberColumns] ; 

3095  memcpy(oldLower,lower,numberColumns*sizeof(double)) ; 

3096  memcpy(oldUpper,upper,numberColumns*sizeof(double)) ; 

3097  #endif 

3098  

3099  double integerTolerance = getDblParam(CbcIntegerTolerance) ; 

3100  for (int i = 0;i<numberColumnCuts;i++) { 

3101  const OsiColCut * thisCut = theseCuts.colCutPtr(i) ; 

3102  const CoinPackedVector & lbs = thisCut>lbs() ; 

3103  const CoinPackedVector & ubs = thisCut>ubs() ; 

3104  int j ; 

3105  int n ; 

3106  const int * which ; 

3107  const double * values ; 

3108  n = lbs.getNumElements() ; 

3109  which = lbs.getIndices() ; 

3110  values = lbs.getElements() ; 

3111  for (j = 0;j<n;j++) { 

3112  int iColumn = which[j] ; 

3113  double value = solution[iColumn] ; 

3114  #if CBC_DEBUG>1 

3115  printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn], 

3116  solution[iColumn],oldUpper[iColumn],values[j]) ; 

3117  #endif 

3118  solver_>setColLower(iColumn,values[j]) ; 

3119  if (value<values[j]integerTolerance) 

3120  violated = 1 ; 

3121  if (values[j]>upper[iColumn]+integerTolerance) { 

3122  // infeasible 

3123  violated = 2 ; 

3124  break ; 

3125  } 

3126  } 

3127  n = ubs.getNumElements() ; 

3128  which = ubs.getIndices() ; 

3129  values = ubs.getElements() ; 

3130  for (j = 0;j<n;j++) { 

3131  int iColumn = which[j] ; 

3132  double value = solution[iColumn] ; 

3133  #if CBC_DEBUG>1 

3134  printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn], 

3135  solution[iColumn],oldUpper[iColumn],values[j]) ; 

3136  #endif 

3137  solver_>setColUpper(iColumn,values[j]) ; 

3138  if (value>values[j]+integerTolerance) 

3139  violated = 1 ; 

3140  if (values[j]<lower[iColumn]integerTolerance) { 

3141  // infeasible 

3142  violated = 2 ; 

3143  break ; 

3144  } 

3145  } 

3146  } 

3147  #ifdef CBC_DEBUG 

3148  delete [] oldLower ; 

3149  delete [] oldUpper ; 

3150  #endif 

3151  } 

3152  /* 

3153  End installation of column cuts. The break here escapes the numberTries 

3154  loop. 

3155  */ 

3156  if (violated == 2!feasible) { 

3157  // infeasible 

3158  feasible = false ; 

3159  violated = 2; 

3160  if (!numberNodes_) 

3161  messageHandler()>message(CBC_INFEAS, 

3162  messages()) 

3163  << CoinMessageEol ; 

3164  break ; 

3165  } 

3166  /* 

3167  Now apply the row (constraint) cuts. This is a bit more work because we need 

3168  to obtain and augment the current basis. 

3169  

3170  TODO: Why do this work, if there are no row cuts? The current basis will do 

3171  just fine. 

3172  */ 

3173  int numberRowsNow = solver_>getNumRows() ; 

3174  assert(numberRowsNow == numberRowsAtStart+lastNumberCuts) ; 

3175  int numberToAdd = theseCuts.sizeRowCuts() ; 

3176  numberNewCuts = lastNumberCuts+numberToAdd ; 

3177  /* 

3178  Get a basis by asking the solver for warm start information. Resize it 

3179  (retaining the basis) so it can accommodate the cuts. 

3180  */ 

3181  delete basis_ ; 

3182  basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

3183  assert(basis_ != NULL); // make sure not volume 

3184  basis_>resize(numberRowsAtStart+numberNewCuts,numberColumns) ; 

3185  /* 

3186  Now actually add the row cuts and reoptimise. 

3187  

3188  Install the cuts in the solver using applyRowCuts and 

3189  augment the basis with the corresponding slack. We also add each row cut to 

3190  the set of row cuts (cuts.insert()) supplied as a parameter. The new basis 

3191  must be set with setWarmStart(). 

3192  

3193  TODO: It's not clear to me why we can't separate this into two sections. 

3194  The first would add the row cuts, and be executed only if row cuts 

3195  need to be installed. The second would call resolve() and would be 

3196  executed if either row or column cuts have been installed. 

3197  

3198  TODO: Seems to me the original code could allocate addCuts with size 0, if 

3199  numberRowCuts was 0 and numberColumnCuts was nonzero. That might 

3200  explain the memory fault noted in the comment by AJK. Unfortunately, 

3201  just commenting out the delete[] results in massive memory leaks. Try 

3202  a revision to separate the row cut case. Why do we need addCuts at 

3203  all? A typing issue, apparently: OsiCut vs. OsiRowCut. 

3204  

3205  TODO: It looks to me as if numberToAdd and numberRowCuts are identical at 

3206  this point. Confirm & get rid of one of them. 

3207  

3208  TODO: Any reason why the three loops can't be consolidated? 

3209  */ 

3210  if (numberRowCuts > 0  numberColumnCuts > 0) 

3211  { if (numberToAdd > 0) 

3212  { int i ; 

3213  // Faster to add all at once 

3214  const OsiRowCut ** addCuts = new const OsiRowCut * [numberToAdd] ; 

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

3216  { addCuts[i] = &theseCuts.rowCut(i) ; } 

3217  solver_>applyRowCuts(numberToAdd,addCuts) ; 

3218  // AJK this caused a memory fault on Win32 

3219  // may be okay now with ** form 

3220  delete [] addCuts ; 

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

3222  { cuts.insert(theseCuts.rowCut(i)) ; } 

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

3224  { basis_>setArtifStatus(numberRowsNow+i, 

3225  CoinWarmStartBasis::basic) ; } 

3226  if (solver_>setWarmStart(basis_) == false) 

3227  { throw CoinError("Fail setWarmStart() after cut installation.", 

3228  "solveWithCuts","CbcModel") ; } } 

3229  feasible = resolve() ; 

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

3231  numberTries=0; // exit 

3232  # ifdef CBC_DEBUG 

3233  printf("Obj value after cuts %g %d rows\n",solver_>getObjValue(), 

3234  solver_>getNumRows()) ; 

3235  if (onOptimalPath && !solver_>isDualObjectiveLimitReached()) 

3236  assert(feasible) ; 

3237  # endif 

3238  } 

3239  /* 

3240  No cuts. Cut short the cut generation (numberTries) loop. 

3241  */ 

3242  else 

3243  { numberTries = 0 ; } 

3244  /* 

3245  If the problem is still feasible, first, call takeOffCuts() to remove cuts 

3246  that are now slack. takeOffCuts() will call the solver to reoptimise if 

3247  that's needed to restore a valid solution. 

3248  

3249  Next, see if we should quit due to diminishing returns: 

3250  * we've tried three rounds of cut generation and we're getting 

3251  insufficient improvement in the objective; or 

3252  * we generated no cuts; or 

3253  * the solver declared optimality with 0 iterations after we added the 

3254  cuts generated in this round. 

3255  If we decide to keep going, prep for the next iteration. 

3256  

3257  It just seems more safe to tell takeOffCuts() to call resolve(), even if 

3258  we're not continuing cut generation. Otherwise code executed between here 

3259  and final disposition of the node will need to be careful not to access the 

3260  lp solution. It can happen that we lose feasibility in takeOffCuts  

3261  numerical jitters when the cutoff bound is epsilon less than the current 

3262  best, and we're evaluating an alternative optimum. 

3263  

3264  TODO: After successive rounds of code motion, there seems no need to 

3265  distinguish between the two checks for aborting the cut generation 

3266  loop. Confirm and clean up. 

3267  */ 

3268  if (feasible) 

3269  { int cutIterations = solver_>getIterationCount() ; 

3270  if (numberOldActiveCuts+numberNewCuts) { 

3271  takeOffCuts(cuts,whichGenerator,numberOldActiveCuts, 

3272  numberNewCuts,resolveAfterTakeOffCuts_) ; 

3273  if (solver_>isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_) 

3274  { feasible = false ; 

3275  # ifdef CBC_DEBUG 

3276  double z = solver_>getObjValue() ; 

3277  double cut = getCutoff() ; 

3278  printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n", 

3279  zcut,z,cut) ; 

3280  # endif 

3281  } 

3282  } 

3283  if (feasible) 

3284  { numberRowsAtStart = numberOldActiveCuts+numberRowsAtContinuous_ ; 

3285  lastNumberCuts = numberNewCuts ; 

3286  if (direction*solver_>getObjValue() < lastObjective+minimumDrop && 

3287  currentPassNumber_ >= 3) 

3288  { numberTries = 0 ; } 

3289  if (numberRowCuts+numberColumnCuts == 0  cutIterations == 0) 

3290  { break ; } 

3291  if (numberTries > 0) 

3292  { reducedCostFix() ; 

3293  lastObjective = direction*solver_>getObjValue() ; 

3294  lower = solver_>getColLower() ; 

3295  upper = solver_>getColUpper() ; 

3296  solution = solver_>getColSolution() ; } } } 

3297  /* 

3298  We've lost feasibility  this node won't be referencing the cuts we've 

3299  been collecting, so decrement the reference counts. 

3300  

3301  TODO: Presumably this is in preparation for backtracking. Seems like it 

3302  should be the `else' off the previous `if'. 

3303  */ 

3304  if (!feasible) 

3305  { int i ; 

3306  for (i = 0;i<currentNumberCuts_;i++) { 

3307  // take off node 

3308  if (addedCuts_[i]) { 

3309  if (!addedCuts_[i]>decrement()) 

3310  delete addedCuts_[i] ; 

3311  addedCuts_[i] = NULL ; 

3312  } 

3313  } 

3314  numberTries = 0 ; 

3315  } 

3316  } while (numberTries>0) ; 

3317  // Reduced cost fix at end 

3318  //reducedCostFix(); 

3319  // If at root node do heuristics 

3320  if (!numberNodes_) { 

3321  // mark so heuristics can tell 

3322  int savePass=currentPassNumber_; 

3323  currentPassNumber_=999999; 

3324  double * newSolution = new double [numberColumns] ; 

3325  double heuristicValue = getCutoff() ; 

3326  int found = 1; // no solution found 

3327  for (int i = 0;i<numberHeuristics_;i++) { 

3328  // see if heuristic will do anything 

3329  double saveValue = heuristicValue ; 

3330  int ifSol = heuristic_[i]>solution(heuristicValue, 

3331  newSolution); 

3332  if (ifSol>0) { 

3333  // better solution found 

3334  found = i ; 

3335  incrementUsed(newSolution); 

3336  } else { 

3337  heuristicValue = saveValue ; 

3338  } 

3339  } 

3340  currentPassNumber_=savePass; 

3341  if (found >= 0) { 

3342  phase_=4; 

3343  incrementUsed(newSolution); 

3344  setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 

3345  } 

3346  delete [] newSolution ; 

3347  } 

3348  // Up change due to cuts 

3349  if (feasible) 

3350  sumChangeObjective2_ += solver_>getObjValue()*solver_>getObjSense() 

3351   objectiveValue ; 

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

3353  //printf("XXb sum obj changed by %g\n",sumChangeObjective2_); 

3354  /* 

3355  End of cut generation loop. 

3356  

3357  Now, consider if we want to disable or adjust the frequency of use for any 

3358  of the cut generators. If the client specified a positive number for 

3359  howOften, it will never change. If the original value was negative, it'll 

3360  be converted to 1000000+howOften, and this value will be adjusted each 

3361  time fullScan is true. Actual cut generation is performed every 

3362  howOften%1000000 nodes; the 1000000 offset is just a convenient way to 

3363  specify that the frequency is adjustable. 

3364  

3365  During cut generation, we recorded the number of cuts produced by each 

3366  generator for this node. For all cuts, whichGenerator records the generator 

3367  that produced a cut. 

3368  

3369  TODO: All this should probably be hidden in a method of the CbcCutGenerator 

3370  class. 

3371  */ 

3372  if (fullScan&&numberCutGenerators_) { 

3373  /* If cuts just at root node then it will probably be faster to 

3374  update matrix and leave all in */ 

3375  bool willBeCutsInTree=false; 

3376  // Root node or every so often  see what to turn off 

3377  int i ; 

3378  double thisObjective = solver_>getObjValue()*direction ; 

3379  double totalCuts = 0.0 ; 

3380  for (i = 0;i<numberCutGenerators_;i++) 

3381  totalCuts += countRowCuts[i] + 5.0*countColumnCuts[i] ; 

3382  if (!numberNodes_) 

3383  handler_>message(CBC_ROOT,messages_) 

3384  <<numberNewCuts 

3385  <<startObjective<<thisObjective 

3386  <<currentPassNumber_ 

3387  <<CoinMessageEol ; 

3388  int * count = new int[numberCutGenerators_] ; 

3389  memset(count,0,numberCutGenerators_*sizeof(int)) ; 

3390  for (i = 0;i<numberNewCuts;i++) { 

3391  int iGenerator = whichGenerator[i]; 

3392  if (iGenerator>=0) 

3393  count[iGenerator]++ ; 

3394  } 

3395  double small = (0.5* totalCuts) / 

3396  ((double) numberCutGenerators_) ; 

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

3398  int howOften = generator_[i]>howOften() ; 

3399  if (howOften<99) 

3400  continue ; 

3401  if (howOften<0howOften >= 1000000) { 

3402  // If small number switch mostly off 

3403  double thisCuts = countRowCuts[i] + 5.0*countColumnCuts[i] ; 

3404  if (!thisCutshowOften == 99) { 

3405  if (howOften == 99) 

3406  howOften = 100 ; 

3407  else 

3408  howOften = 1000000+SCANCUTS; // wait until next time 

3409  } else if (thisCuts<small) { 

3410  int k = (int) sqrt(small/thisCuts) ; 

3411  howOften = k+1000000 ; 

3412  } else { 

3413  howOften = 1+1000000 ; 

3414  } 

3415  // If cuts useless switch off 

3416  if (numberNodes_>=10&&sumChangeObjective1_>1.0e2*(sumChangeObjective2_+1.0e12)) { 

3417  howOften = 1000000+SCANCUTS; // wait until next time 

3418  //printf("switch off cut %d due to lack of use\n",i); 

3419  } 

3420  } 

3421  if (howOften>=0&&generator_[i]>generator()>mayGenerateRowCutsInTree()) 

3422  willBeCutsInTree=true; 

3423  

3424  generator_[i]>setHowOften(howOften) ; 

3425  if (howOften>=1000000&&howOften<2000000&&0) { 

3426  // Go to depth 

3427  int bias=1; 

3428  if (howOften==1+1000000) 

3429  generator_[i]>setWhatDepth(bias+1); 

3430  else if (howOften<=10+1000000) 

3431  generator_[i]>setWhatDepth(bias+2); 

3432  else 

3433  generator_[i]>setWhatDepth(bias+1000); 

3434  } 

3435  int newFrequency = generator_[i]>howOften()%1000000 ; 

3436  // increment cut counts 

3437  generator_[i]>incrementNumberCutsInTotal(countRowCuts[i]); 

3438  generator_[i]>incrementNumberCutsActive(count[i]); 

3439  if (handler_>logLevel()>1!numberNodes_) { 

3440  handler_>message(CBC_GENERATOR,messages_) 

3441  <<i 

3442  <<generator_[i]>cutGeneratorName() 

3443  <<countRowCuts[i]<<count[i] 

3444  <<countColumnCuts[i]; 

3445  handler_>printing(!numberNodes_&&generator_[i]>timing()) 

3446  <<generator_[i]>timeInCutGenerator(); 

3447  handler_>message() 

3448  <<newFrequency 

3449  <<CoinMessageEol ; 

3450  } 

3451  } 

3452  delete [] count ; 

3453  if( !numberNodes_) { 

3454  if( !willBeCutsInTree) { 

3455  // Take off cuts 

3456  cuts = OsiCuts(); 

3457  numberNewCuts=0; 

3458  // update size of problem 

3459  numberRowsAtContinuous_ = solver_>getNumRows() ; 

3460  #ifdef COIN_USE_CLP 

3461  OsiClpSolverInterface * clpSolver 

3462  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

3463  if (clpSolver) { 

3464  // Maybe solver might like to know only column bounds will change 

3465  //int options = clpSolver>specialOptions(); 

3466  //clpSolver>setSpecialOptions(options128); 

3467  clpSolver>synchronizeModel(); 

3468  } 

3469  #endif 

3470  } else { 

3471  #ifdef COIN_USE_CLP 

3472  OsiClpSolverInterface * clpSolver 

3473  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

3474  if (clpSolver) { 

3475  // make sure factorization can't carry over 

3476  int options = clpSolver>specialOptions(); 

3477  clpSolver>setSpecialOptions(options&(~8)); 

3478  } 

3479  #endif 

3480  } 

3481  } 

3482  } else if (numberCutGenerators_) { 

3483  int i; 

3484  // add to counts anyway 

3485  for (i = 0;i<numberCutGenerators_;i++) 

3486  generator_[i]>incrementNumberCutsInTotal(countRowCuts[i]); 

3487  // What if not feasible as cuts may have helped 

3488  if (feasible) { 

3489  for (i = 0;i<numberNewCuts;i++) { 

3490  int iGenerator = whichGenerator[i]; 

3491  if (iGenerator>=0) 

3492  generator_[iGenerator]>incrementNumberCutsActive(); 

3493  } 

3494  } 

3495  } 

3496  

3497  delete [] countRowCuts ; 

3498  delete [] countColumnCuts ; 

3499  

3500  

3501  #ifdef CHECK_CUT_COUNTS 

3502  if (feasible) 

3503  { delete basis_ ; 

3504  basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

3505  printf("solveWithCuts: Number of rows at end (only active cuts) %d\n", 

3506  numberRowsAtContinuous_+numberNewCuts+numberOldActiveCuts) ; 

3507  basis_>print() ; } 

3508  #endif 

3509  #ifdef CBC_DEBUG 

3510  if (onOptimalPath && !solver_>isDualObjectiveLimitReached()) 

3511  assert(feasible) ; 

3512  #endif 

3513  

3514  return feasible ; } 

3515  

3516  

3517  /* 

3518  Remove slack cuts. We obtain a basis and scan it. Cuts with basic slacks 

3519  are purged. If any cuts are purged, resolve() is called to restore the 

3520  solution held in the solver. If resolve() pivots, there's the possibility 

3521  that a slack may be pivoted in (trust me :), so the process iterates. 

3522  Setting allowResolve to false will suppress reoptimisation (but see note 

3523  below). 

3524  

3525  At the level of the solver's constraint system, loose cuts are really 

3526  deleted. There's an implicit assumption that deleteRows will also update 

3527  the active basis in the solver. 

3528  

3529  At the level of nodes and models, it's more complicated. 

3530  

3531  New cuts exist only in the collection of cuts passed as a parameter. They 

3532  are deleted from the collection and that's the end of them. 

3533  

3534  Older cuts have made it into addedCuts_. Two separate actions are needed. 

3535  The reference count for the CbcCountRowCut object is decremented. If this 

3536  count falls to 0, the node which owns the cut is located, the reference to 

3537  the cut is removed, and then the cut object is destroyed (courtesy of the 

3538  CbcCountRowCut destructor). We also need to set the addedCuts_ entry to 

3539  NULL. This is important so that when it comes time to generate basis edits 

3540  we can tell this cut was dropped from the basis during processing of the 

3541  node. 

3542  

3543  NOTE: In general, it's necessary to call resolve() after purging slack 

3544  cuts. Deleting constraints constitutes a change in the problem, and 

3545  an OSI is not required to maintain a valid solution when the problem 

3546  is changed. But ... it's really useful to maintain the active basis, 

3547  and the OSI is supposed to do that. (Yes, it's splitting hairs.) In 

3548  some places, it's possible to know that the solution will never be 

3549  consulted after this call, only the basis. (E.g., this routine is 

3550  called as a last act before generating info to place the node in the 

3551  live set.) For such use, set allowResolve to false. 

3552  

3553  TODO: No real harm would be done if we just ignored the rare occasion when 

3554  the call to resolve() pivoted a slack back into the basis. It's a 

3555  minor inefficiency, at worst. But it does break assertions which 

3556  check that there are no loose cuts in the basis. It might be better 

3557  to remove the assertions. 

3558  */ 

3559  

3560  void 

3561  CbcModel::takeOffCuts (OsiCuts &newCuts, int *whichGenerator, 

3562  int &numberOldActiveCuts, int &numberNewCuts, 

3563  bool allowResolve) 

3564  

3565  { // int resolveIterations = 0 ; 

3566  int firstOldCut = numberRowsAtContinuous_ ; 

3567  int totalNumberCuts = numberNewCuts+numberOldActiveCuts ; 

3568  int *solverCutIndices = new int[totalNumberCuts] ; 

3569  int *newCutIndices = new int[numberNewCuts] ; 

3570  const CoinWarmStartBasis* ws ; 

3571  CoinWarmStartBasis::Status status ; 

3572  bool needPurge = true ; 

3573  /* 

3574  The outer loop allows repetition of purge in the event that reoptimisation 

3575  changes the basis. To start an iteration, clear the deletion counts and grab 

3576  the current basis. 

3577  */ 

3578  while (needPurge) 

3579  { int numberNewToDelete = 0 ; 

3580  int numberOldToDelete = 0 ; 

3581  int i ; 

3582  ws = dynamic_cast<const CoinWarmStartBasis*>(solver_>getWarmStart()) ; 

3583  /* 

3584  Scan the basis entries of the old cuts generated prior to this round of cut 

3585  generation. Loose cuts are `removed' by decrementing their reference count 

3586  and setting the addedCuts_ entry to NULL. (If the reference count falls to 

3587  0, they're really deleted. See CbcModel and CbcCountRowCut doc'n for 

3588  principles of cut handling.) 

3589  */ 

3590  int oldCutIndex = 0 ; 

3591  for (i = 0 ; i < numberOldActiveCuts ; i++) 

3592  { status = ws>getArtifStatus(i+firstOldCut) ; 

3593  while (!addedCuts_[oldCutIndex]) oldCutIndex++ ; 

3594  assert(oldCutIndex < currentNumberCuts_) ; 

3595  // always leave if from nextRowCut_ 

3596  if (status == CoinWarmStartBasis::basic&& 

3597  addedCuts_[oldCutIndex]>effectiveness()!=COIN_DBL_MAX) 

3598  { solverCutIndices[numberOldToDelete++] = i+firstOldCut ; 

3599  if (addedCuts_[oldCutIndex]>decrement() == 0) 

3600  delete addedCuts_[oldCutIndex] ; 

3601  addedCuts_[oldCutIndex] = NULL ; 

3602  oldCutIndex++ ; } 

3603  else 

3604  { oldCutIndex++ ; } } 

3605  /* 

3606  Scan the basis entries of the new cuts generated with this round of cut 

3607  generation. At this point, newCuts is the only record of the new cuts, so 

3608  when we delete loose cuts from newCuts, they're really gone. newCuts is a 

3609  vector, so it's most efficient to compress it (eraseRowCut) from back to 

3610  front. 

3611  */ 

3612  int firstNewCut = firstOldCut+numberOldActiveCuts ; 

3613  int k = 0 ; 

3614  for (i = 0 ; i < numberNewCuts ; i++) 

3615  { status = ws>getArtifStatus(i+firstNewCut) ; 

3616  if (status == CoinWarmStartBasis::basic&&whichGenerator[i]!=2) 

3617  { solverCutIndices[numberNewToDelete+numberOldToDelete] = i+firstNewCut ; 

3618  newCutIndices[numberNewToDelete++] = i ; } 

3619  else 

3620  { // save which generator did it 

3621  whichGenerator[k++] = whichGenerator[i] ; } } 

3622  delete ws ; 

3623  for (i = numberNewToDelete1 ; i >= 0 ; i) 

3624  { int iCut = newCutIndices[i] ; 

3625  newCuts.eraseRowCut(iCut) ; } 

3626  /* 

3627  Did we delete anything? If so, delete the cuts from the constraint system 

3628  held in the solver and reoptimise unless we're forbidden to do so. If the 

3629  call to resolve() results in pivots, there's the possibility we again have 

3630  basic slacks. Repeat the purging loop. 

3631  */ 

3632  if (numberNewToDelete+numberOldToDelete > 0) 

3633  { solver_>deleteRows(numberNewToDelete+numberOldToDelete, 

3634  solverCutIndices) ; 

3635  numberNewCuts = numberNewToDelete ; 

3636  numberOldActiveCuts = numberOldToDelete ; 

3637  # ifdef CBC_DEBUG 

3638  printf("takeOffCuts: purged %d+%d cuts\n", numberOldToDelete, 

3639  numberNewToDelete ); 

3640  # endif 

3641  if (allowResolve) 

3642  { 

3643  phase_=3; 

3644  // can do quick optimality check 

3645  int easy=2; 

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

3647  solver_>resolve() ; 

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

3649  if (solver_>getIterationCount() == 0) 

3650  { needPurge = false ; } 

3651  # ifdef CBC_DEBUG 

3652  else 

3653  { printf( "Repeating purging loop. %d iters.\n", 

3654  solver_>getIterationCount()); 

3655  # endif 

3656  } 

3657  else 

3658  { needPurge = false ; } } 

3659  else 

3660  { needPurge = false ; } } 

3661  /* 

3662  Clean up and return. 

3663  */ 

3664  delete [] solverCutIndices ; 

3665  delete [] newCutIndices ; 

3666  } 

3667  

3668  

3669  

3670  bool 

3671  CbcModel::resolve() 

3672  { 

3673  // We may have deliberately added in violated cuts  check to avoid message 

3674  int iRow; 

3675  int numberRows = solver_>getNumRows(); 

3676  const double * rowLower = solver_>getRowLower(); 

3677  const double * rowUpper = solver_>getRowUpper(); 

3678  bool feasible=true; 

3679  for (iRow= numberRowsAtContinuous_;iRow<numberRows;iRow++) { 

3680  if (rowLower[iRow]>rowUpper[iRow]+1.0e8) 

3681  feasible=false; 

3682  } 

3683  // Can't happen if strong branching as would have been found before 

3684  if (!numberStrong_&&numberObjects_>numberIntegers_) { 

3685  int iColumn; 

3686  int numberColumns = solver_>getNumCols(); 

3687  const double * columnLower = solver_>getColLower(); 

3688  const double * columnUpper = solver_>getColUpper(); 

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

3690  if (columnLower[iColumn]>columnUpper[iColumn]+1.0e5) 

3691  feasible=false; 

3692  } 

3693  } 

3694  /* 

3695  Reoptimize. Consider the possibility that we should fathom on bounds. But be 

3696  careful  where the objective takes on integral values, we may want to keep 

3697  a solution where the objective is right on the cutoff. 

3698  */ 

3699  if (feasible) 

3700  { 

3701  solver_>resolve() ; 

3702  numberIterations_ += solver_>getIterationCount() ; 

3703  feasible = (solver_>isProvenOptimal() && 

3704  !solver_>isDualObjectiveLimitReached()) ; } 

3705  if (!feasible&& continuousObjective_ <1.0e30) { 

3706  // at root node  double double check 

3707  bool saveTakeHint; 

3708  OsiHintStrength saveStrength; 

3709  solver_>getHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength); 

3710  if (saveTakeHintsaveStrength==OsiHintIgnore) { 

3711  solver_>setHintParam(OsiDoDualInResolve,false,OsiHintDo) ; 

3712  solver_>resolve(); 

3713  solver_>setHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength); 

3714  numberIterations_ += solver_>getIterationCount() ; 

3715  feasible = solver_>isProvenOptimal(); 

3716  } 

3717  } 

3718  return feasible ; } 

3719  

3720  

3721  /* Set up objects. Only do ones whose length is in range. 

3722  If makeEquality true then a new model may be returned if 

3723  modifications had to be made, otherwise "this" is returned. 

3724  

3725  Could use Probing at continuous to extend objects 

3726  */ 

3727  CbcModel * 

3728  CbcModel::findCliques(bool makeEquality, 

3729  int atLeastThisMany, int lessThanThis, int defaultValue) 

3730  { 

3731  // No objects are allowed to exist 

3732  assert(numberObjects_==numberIntegers_!numberObjects_); 

3733  CoinPackedMatrix matrixByRow(*solver_>getMatrixByRow()); 

3734  int numberRows = solver_>getNumRows(); 

3735  int numberColumns = solver_>getNumCols(); 

3736  

3737  // We may want to add columns 

3738  int numberSlacks=0; 

3739  int * rows = new int[numberRows]; 

3740  double * element =new double[numberRows]; 

3741  

3742  int iRow; 

3743  

3744  findIntegers(true); 

3745  numberObjects_=numberIntegers_; 

3746  

3747  int numberCliques=0; 

3748  CbcObject ** object = new CbcObject * [numberRows]; 

3749  int * which = new int[numberIntegers_]; 

3750  char * type = new char[numberIntegers_]; 

3751  int * lookup = new int[numberColumns]; 

3752  int i; 

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

3754  lookup[i]=1; 

3755  for (i=0;i<numberIntegers_;i++) 

3756  lookup[integerVariable_[i]]=i; 

3757  

3758  // Statistics 

3759  int totalP1=0,totalM1=0; 

3760  int numberBig=0,totalBig=0; 

3761  int numberFixed=0; 

3762  

3763  // Row copy 

3764  const double * elementByRow = matrixByRow.getElements(); 

3765  const int * column = matrixByRow.getIndices(); 

3766  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts(); 

3767  const int * rowLength = matrixByRow.getVectorLengths(); 

3768  

3769  // Column lengths for slacks 

3770  const int * columnLength = solver_>getMatrixByCol()>getVectorLengths(); 

3771  

3772  const double * lower = getColLower(); 

3773  const double * upper = getColUpper(); 

3774  const double * rowLower = getRowLower(); 

3775  const double * rowUpper = getRowUpper(); 

3776  

3777  for (iRow=0;iRow<numberRows;iRow++) { 

3778  int numberP1=0, numberM1=0; 

3779  int j; 

3780  double upperValue=rowUpper[iRow]; 

3781  double lowerValue=rowLower[iRow]; 

3782  bool good=true; 

3783  int slack = 1; 

3784  for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) { 

3785  int iColumn = column[j]; 

3786  int iInteger=lookup[iColumn]; 

3787  if (upper[iColumn]lower[iColumn]<1.0e8) { 

3788  // fixed 

3789  upperValue = lower[iColumn]*elementByRow[j]; 

3790  lowerValue = lower[iColumn]*elementByRow[j]; 

3791  continue; 

3792  } else if (upper[iColumn]!=1.0lower[iColumn]!=0.0) { 

3793  good = false; 

3794  break; 

3795  } else if (iInteger<0) { 

3796  good = false; 

3797  break; 

3798  } else { 

3799  if (columnLength[iColumn]==1) 

3800  slack = iInteger; 

3801  } 

3802  if (fabs(elementByRow[j])!=1.0) { 

3803  good=false; 

3804  break; 

3805  } else if (elementByRow[j]>0.0) { 

3806  which[numberP1++]=iInteger; 

3807  } else { 

3808  numberM1++; 

3809  which[numberIntegers_numberM1]=iInteger; 

3810  } 

3811  } 

3812  int iUpper = (int) floor(upperValue+1.0e5); 

3813  int iLower = (int) ceil(lowerValue1.0e5); 

3814  int state=0; 

3815  if (upperValue<1.0e6) { 

3816  if (iUpper==1numberM1) 

3817  state=1; 

3818  else if (iUpper==numberM1) 

3819  state=2; 

3820  else if (iUpper<numberM1) 

3821  state=3; 

3822  } 

3823  if (!state&&lowerValue>1.0e6) { 

3824  if (iLower==1numberP1) 

3825  state=1; 

3826  else if (iLower==numberP1) 

3827  state=2; 

3828  else if (iLower<numberP1) 

3829  state=3; 

3830  } 

3831  if (good&&state) { 

3832  if (abs(state)==3) { 

3833  // infeasible 

3834  numberObjects_=1; 

3835  break; 

3836  } else if (abs(state)==2) { 

3837  // we can fix all 

3838  numberFixed += numberP1+numberM1; 

3839  if (state>0) { 

3840  // fix all +1 at 0, 1 at 1 

3841  for (i=0;i<numberP1;i++) 

3842  solver_>setColUpper(integerVariable_[which[i]],0.0); 

3843  for (i=0;i<numberM1;i++) 

3844  solver_>setColLower(integerVariable_[which[numberIntegers_i1]], 

3845  1.0); 

3846  } else { 

3847  // fix all +1 at 1, 1 at 0 

3848  for (i=0;i<numberP1;i++) 

3849  solver_>setColLower(integerVariable_[which[i]],1.0); 

3850  for (i=0;i<numberM1;i++) 

3851  solver_>setColUpper(integerVariable_[which[numberIntegers_i1]], 

3852  0.0); 

3853  } 

3854  } else { 

3855  int length = numberP1+numberM1; 

3856  if (length >= atLeastThisMany&&length<lessThanThis) { 

3857  // create object 

3858  bool addOne=false; 

3859  int objectType; 

3860  if (iLower==iUpper) { 

3861  objectType=1; 

3862  } else { 

3863  if (makeEquality) { 

3864  objectType=1; 

3865  element[numberSlacks]=state; 

3866  rows[numberSlacks++]=iRow; 

3867  addOne=true; 

3868  } else { 

3869  objectType=0; 

3870  } 

3871  } 

3872  if (state>0) { 

3873  totalP1 += numberP1; 

3874  totalM1 += numberM1; 

3875  for (i=0;i<numberP1;i++) 

3876  type[i]=1; 

3877  for (i=0;i<numberM1;i++) { 

3878  which[numberP1]=which[numberIntegers_i1]; 

3879  type[numberP1++]=0; 

3880  } 

3881  } else { 

3882  totalP1 += numberM1; 

3883  totalM1 += numberP1; 

3884  for (i=0;i<numberP1;i++) 

3885  type[i]=0; 

3886  for (i=0;i<numberM1;i++) { 

3887  which[numberP1]=which[numberIntegers_i1]; 

3888  type[numberP1++]=1; 

3889  } 

3890  } 

3891  if (addOne) { 

3892  // add in slack 

3893  which[numberP1]=numberIntegers_+numberSlacks1; 

3894  slack = numberP1; 

3895  type[numberP1++]=1; 

3896  } else if (slack >= 0) { 

3897  for (i=0;i<numberP1;i++) { 

3898  if (which[i]==slack) { 

3899  slack=i; 

3900  } 

3901  } 

3902  } 

3903  object[numberCliques] = new CbcClique(this,objectType,numberP1, 

3904  which,type, 

3905  1000000+numberCliques,slack); 

3906  numberCliques++; 

3907  } else if (numberP1+numberM1 >= lessThanThis) { 

3908  // too big 

3909  numberBig++; 

3910  totalBig += numberP1+numberM1; 

3911  } 

3912  } 

3913  } 

3914  } 

3915  delete [] which; 

3916  delete [] type; 

3917  delete [] lookup; 

3918  if (numberCliques<0) { 

3919  printf("*** Problem infeasible\n"); 

3920  } else { 

3921  if (numberCliques) 

3922  printf("%d cliques of average size %g found, %d P1, %d M1\n", 

3923  numberCliques, 

3924  ((double)(totalP1+totalM1))/((double) numberCliques), 

3925  totalP1,totalM1); 

3926  else 

3927  printf("No cliques found\n"); 

3928  if (numberBig) 

3929  printf("%d large cliques ( >= %d) found, total %d\n", 

3930  numberBig,lessThanThis,totalBig); 

3931  if (numberFixed) 

3932  printf("%d variables fixed\n",numberFixed); 

3933  } 

3934  if (numberCliques>0&&numberSlacks&&makeEquality) { 

3935  printf("adding %d integer slacks\n",numberSlacks); 

3936  // add variables to make equality rows 

3937  int * temp = new int[numberIntegers_+numberSlacks]; 

3938  memcpy(temp,integerVariable_,numberIntegers_*sizeof(int)); 

3939  // Get new model 

3940  CbcModel * newModel = new CbcModel(*this); 

3941  OsiSolverInterface * newSolver = newModel>solver(); 

3942  for (i=0;i<numberSlacks;i++) { 

3943  temp[i+numberIntegers_]=i+numberColumns; 

3944  int iRow = rows[i]; 

3945  double value = element[i]; 

3946  double lowerValue = 0.0; 

3947  double upperValue = 1.0; 

3948  double objValue = 0.0; 

3949  CoinPackedVector column(1,&iRow,&value); 

3950  newSolver>addCol(column,lowerValue,upperValue,objValue); 

3951  // set integer 

3952  newSolver>setInteger(numberColumns+i); 

3953  if (value >0) 

3954  newSolver>setRowLower(iRow,rowUpper[iRow]); 

3955  else 

3956  newSolver>setRowUpper(iRow,rowLower[iRow]); 

3957  } 

3958  // replace list of integers 

3959  for (i=0;i<newModel>numberObjects_;i++) 

3960  delete newModel>object_[i]; 

3961  newModel>numberObjects_ = 0; 

3962  delete [] newModel>object_; 

3963  newModel>object_=NULL; 

3964  newModel>findIntegers(true); //Set up all integer objects 

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

3966  newModel>modifiableObject(i)>setPriority(object_[i]>priority()); 

3967  } 

3968  if (originalColumns_) { 

3969  // old model had originalColumns 

3970  delete [] newModel>originalColumns_; 

3971  newModel>originalColumns_ = new int[numberColumns+numberSlacks]; 

3972  memcpy(newModel>originalColumns_,originalColumns_,numberColumns*sizeof(int)); 

3973  // mark as not in previous model 

3974  for (i=numberColumns;i<numberColumns+numberSlacks;i++) 

3975  newModel>originalColumns_[i]=1; 

3976  } 

3977  delete [] rows; 

3978  delete [] element; 

3979  newModel>addObjects(numberCliques,object); 

3980  for (;i<numberCliques;i++) 

3981  delete object[i]; 

3982  delete [] object; 

3983  newModel>synchronizeModel(); 

3984  return newModel; 

3985  } else { 

3986  if (numberCliques>0) { 

3987  addObjects(numberCliques,object); 

3988  for (;i<numberCliques;i++) 

3989  delete object[i]; 

3990  synchronizeModel(); 

3991  } 

3992  delete [] object; 

3993  delete [] rows; 

3994  delete [] element; 

3995  return this; 

3996  } 

3997  } 

3998  

3999  /* 

4000  Set branching priorities. 

4001  

4002  Setting integer priorities looks pretty robust; the call to findIntegers 

4003  makes sure that SimpleInteger objects are in place. Setting priorities for 

4004  other objects is entirely dependent on their existence, and the routine may 

4005  quietly fail in several directions. 

4006  */ 

4007  

4008  void 

4009  CbcModel::passInPriorities (const int * priorities, 

4010  bool ifObject) 

4011  { 

4012  findIntegers(false); 

4013  int i; 

4014  if (priorities) { 

4015  int i0=0; 

4016  int i1=numberObjects_1; 

4017  if (ifObject) { 

4018  for (i=numberIntegers_;i<numberObjects_;i++) { 

4019  object_[i]>setPriority(priorities[inumberIntegers_]); 

4020  } 

4021  i0=numberIntegers_; 

4022  } else { 

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

4024  object_[i]>setPriority(priorities[i]); 

4025  } 

4026  i1=numberIntegers_1; 

4027  } 

4028  messageHandler()>message(CBC_PRIORITY, 

4029  messages()) 

4030  << i0<<i1<<numberObjects_ << CoinMessageEol ; 

4031  } 

4032  } 

4033  

4034  // Delete all object information 

4035  void 

4036  CbcModel::deleteObjects() 

4037  { 

4038  int i; 

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

4040  delete object_[i]; 

4041  delete [] object_; 

4042  object_ = NULL; 

4043  numberObjects_=0; 

4044  findIntegers(true); 

4045  } 

4046  

4047  /*! 

4048  Ensure all attached objects (CbcObjects, heuristics, and cut 

4049  generators) point to this model. 

4050  */ 

4051  void CbcModel::synchronizeModel() 

4052  { 

4053  int i; 

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

4055  heuristic_[i]>setModel(this); 

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

4057  object_[i]>setModel(this); 

4058  for (i=0;i<numberCutGenerators_;i++) 

4059  generator_[i]>refreshModel(this); 

4060  } 

4061  

4062  // Fill in integers and create objects 

4063  

4064  /** 

4065  The routine first does a scan to count the number of integer variables. 

4066  It then creates an array, integerVariable_, to store the indices of the 

4067  integer variables, and an array of `objects', one for each variable. 

4068  

4069  The scan is repeated, this time recording the index of each integer 

4070  variable in integerVariable_, and creating an CbcSimpleInteger object that 

4071  contains information about the integer variable. Initially, this is just 

4072  the index and upper & lower bounds. 

4073  

4074  \todo 

4075  Note the assumption in cbc that the first numberIntegers_ objects are 

4076  CbcSimpleInteger. In particular, the code which handles the startAgain 

4077  case assumes that if the object_ array exists it can simply replace the first 

4078  numberInteger_ objects. This is arguably unsafe. 

4079  

4080  I am going to reorder if necessary 

4081  */ 

4082  

4083  void 

4084  CbcModel::findIntegers(bool startAgain) 

4085  { 

4086  assert(solver_); 

4087  /* 

4088  No need to do this if we have previous information, unless forced to start 

4089  over. 

4090  */ 

4091  if (numberIntegers_&&!startAgain&&object_) 

4092  return; 

4093  /* 

4094  Clear out the old integer variable list, then count the number of integer 

4095  variables. 

4096  */ 

4097  delete [] integerVariable_; 

4098  numberIntegers_=0; 

4099  int numberColumns = getNumCols(); 

4100  int iColumn; 

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

4102  if (isInteger(iColumn)) 

4103  numberIntegers_++; 

4104  } 

4105  // Find out how many old noninteger objects there are 

4106  int nObjects=0; 

4107  CbcObject ** oldObject = object_; 

4108  int iObject; 

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

4110  CbcSimpleInteger * obj = 

4111  dynamic_cast <CbcSimpleInteger *>(oldObject[iObject]) ; 

4112  if (obj) 

4113  delete oldObject[iObject]; 

4114  else 

4115  oldObject[nObjects++]=oldObject[iObject]; 

4116  } 

4117  

4118  /* 

4119  Found any? Allocate an array to hold the indices of the integer variables. 

4120  Make a large enough array for all objects 

4121  */ 

4122  object_ = new CbcObject * [numberIntegers_+nObjects]; 

4123  numberObjects_=numberIntegers_+nObjects;; 

4124  integerVariable_ = new int [numberIntegers_]; 

4125  /* 

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

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

4128  lower bounds. 

4129  */ 

4130  numberIntegers_=0; 

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

4132  if(isInteger(iColumn)) { 

4133  object_[numberIntegers_] = 

4134  new CbcSimpleInteger(this,numberIntegers_,iColumn); 

4135  integerVariable_[numberIntegers_++]=iColumn; 

4136  } 

4137  } 

4138  // Now append other objects 

4139  memcpy(object_+numberIntegers_,oldObject,nObjects*sizeof(CbcObject *)); 

4140  // Delete old array (just array) 

4141  delete [] oldObject; 

4142  

4143  if (!numberObjects_) 

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

4145  } 

4146  /* If numberBeforeTrust >0 then we are going to use CbcBranchDynamic. 

4147  Scan and convert CbcSimpleInteger objects 

4148  */ 

4149  void 

4150  CbcModel::convertToDynamic() 

4151  { 

4152  int iObject; 

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

4154  CbcSimpleInteger * obj1 = 

4155  dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ; 

4156  CbcSimpleIntegerDynamicPseudoCost * obj2 = 

4157  dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[iObject]) ; 

4158  if (obj1&&!obj2) { 

4159  // replace 

4160  int iColumn = obj1>columnNumber(); 

4161  delete object_[iObject]; 

4162  CbcSimpleIntegerDynamicPseudoCost * newObject = 

4163  new CbcSimpleIntegerDynamicPseudoCost(this,iObject,iColumn,0.3); 

4164  newObject>setNumberBeforeTrust(numberBeforeTrust_); 

4165  object_[iObject] = newObject; 

4166  } 

4167  } 

4168  if (branchingMethod_&&(branchingMethod_>whichMethod()&1)==0) { 

4169  // Need a method which can do better 

4170  branchingMethod_=NULL; 

4171  } 

4172  } 

4173  

4174  /* Add in any object information (objects are cloned  owner can delete 

4175  originals */ 

4176  void 

4177  CbcModel::addObjects(int numberObjects, CbcObject ** objects) 

4178  { 

4179  // If integers but not enough objects fudge 

4180  if (numberIntegers_>numberObjects_) 

4181  findIntegers(true); 

4182  /* But if incoming objects inherit from simple integer we just want 

4183  to replace */ 

4184  int numberColumns = solver_>getNumCols(); 

4185  /** mark is 1 if not integer, >=0 if using existing simple integer and 

4186  >=numberColumns if using new integer */ 

4187  int * mark = new int[numberColumns]; 

4188  int i; 

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

4190  mark[i]=1; 

4191  int newNumberObjects = numberObjects; 

4192  int newIntegers=0; 

4193  for (i=0;i<numberObjects;i++) { 

4194  CbcSimpleInteger * obj = 

4195  dynamic_cast <CbcSimpleInteger *>(objects[i]) ; 

4196  if (obj) { 

4197  int iColumn = obj>columnNumber(); 

4198  mark[iColumn]=i+numberColumns; 

4199  newIntegers++; 

4200  } 

4201  } 

4202  // and existing 

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

4204  CbcSimpleInteger * obj = 

4205  dynamic_cast <CbcSimpleInteger *>(object_[i]) ; 

4206  if (obj) { 

4207  int iColumn = obj>columnNumber(); 

4208  if (mark[iColumn]<0) { 

4209  newIntegers++; 

4210  newNumberObjects++; 

4211  mark[iColumn]=i; 

4212  } 

4213  } 

4214  } 

4215  delete [] integerVariable_; 

4216  integerVariable_=NULL; 

4217  if (newIntegers!=numberIntegers_) 

4218  printf("changing number of integers from %d to %d\n", 

4219  numberIntegers_,newIntegers); 

4220  numberIntegers_ = newIntegers; 

4221  integerVariable_ = new int [numberIntegers_]; 

4222  CbcObject ** temp = new CbcObject * [newNumberObjects]; 

4223  // Put integers first 

4224  newIntegers=0; 

4225  numberIntegers_=0; 

4226  for (i=0;i<numberColumns;i++) { 

4227  int which = mark[i]; 

4228  if (which>=0) { 

4229  if (!isInteger(i)) { 

4230  newIntegers++; 

4231  solver_>setInteger(i); 

4232  } 

4233  if (which<numberColumns) { 

4234  temp[numberIntegers_]=object_[which]; 

4235  object_[which]=NULL; 

4236  } else { 

4237  temp[numberIntegers_]=objects[whichnumberColumns]>clone(); 

4238  temp[numberIntegers_]>setModel(this); 

4239  } 

4240  integerVariable_[numberIntegers_++]=i; 

4241  } 

4242  } 

4243  if (newIntegers) 

4244  printf("%d variables were declared integer\n",newIntegers); 

4245  int n=numberIntegers_; 

4246  // Now rest of old 

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

4248  if (object_[i]) { 

4249  CbcSimpleInteger * obj = 

4250  dynamic_cast <CbcSimpleInteger *>(object_[i]) ; 

4251  if (obj) { 

4252  delete object_[i]; 

4253  } else { 

4254  temp[n++]=object_[i]; 

4255  } 

4256  } 

4257  } 

4258  // and rest of new 

4259  for (i=0;i<numberObjects;i++) { 

4260  CbcSimpleInteger * obj = 

4261  dynamic_cast <CbcSimpleInteger *>(objects[i]) ; 

4262  if (!obj) { 

4263  temp[n]=objects[i]>clone(); 

4264  temp[n++]>setModel(this); 

4265  } 

4266  } 

4267  delete [] mark; 

4268  delete [] object_; 

4269  object_ = temp; 

4270  assert (n==newNumberObjects); 

4271  numberObjects_ = newNumberObjects; 

4272  } 

4273  

4274  /** 

4275  This routine sets the objective cutoff value used for fathoming and 

4276  determining monotonic variables. 

4277  

4278  If the fathoming discipline is strict, a small tolerance is added to the 

4279  new cutoff. This avoids problems due to roundoff when the target value 

4280  is exact. The common example would be an IP with only integer variables in 

4281  the objective. If the target is set to the exact value z of the optimum, 

4282  it's possible to end up fathoming an ancestor of the solution because the 

4283  solver returns z+epsilon. 

4284  

4285  Determining if strict fathoming is needed is best done by analysis. 

4286  In cbc, that's analyseObjective. The default is false. 

4287  

4288  In cbc we always minimize so add epsilon 

4289  */ 

4290  

4291  void CbcModel::setCutoff (double value) 

4292  

4293  { 

4294  #if 0 

4295  double tol = 0 ; 

4296  int fathomStrict = getIntParam(CbcFathomDiscipline) ; 

4297  if (fathomStrict == 1) 

4298  { solver_>getDblParam(OsiDualTolerance,tol) ; 

4299  tol = tol*(1+fabs(value)) ; 

4300  

4301  value += tol ; } 

4302  #endif 

4303  // Solvers know about direction 

4304  double direction = solver_>getObjSense(); 

4305  solver_>setDblParam(OsiDualObjectiveLimit,value*direction); } 

4306  

4307  

4308  

4309  /* 

4310  Call this to really test if a valid solution can be feasible. The cutoff is 

4311  passed in as a parameter so that we don't need to worry here after swapping 

4312  solvers. The solution is assumed to be numberColumns in size. If 

4313  fixVariables is true then the bounds of the continuous solver are updated. 

4314  The routine returns the objective value determined by reoptimizing from 

4315  scratch. If the solution is rejected, this will be worse than the cutoff. 

4316  

4317  TODO: There's an issue with getting the correct cutoff value: We update the 

4318  cutoff in the regular solver, but not in continuousSolver_. But our only 

4319  use for continuousSolver_ is verifying candidate solutions. Would it 

4320  make sense to update the cutoff? Then we wouldn't need to step around 

4321  isDualObjectiveLimitReached(). 

4322  */ 

4323  double 

4324  CbcModel::checkSolution (double cutoff, const double *solution, 

4325  bool fixVariables) 

4326  

4327  { int numberColumns = solver_>getNumCols(); 

4328  

4329  /* 

4330  Grab the continuous solver (the pristine copy of the problem, made before 

4331  starting to work on the root node). Save the bounds on the variables. 

4332  Install the solution passed as a parameter, and copy it to the model's 

4333  currentSolution_. 

4334  

4335  TODO: This is a beltandsuspenders approach. Once the code has settled 

4336  a bit, we can cast a critical eye here. 

4337  */ 

4338  OsiSolverInterface * saveSolver = solver_; 

4339  if (continuousSolver_) 

4340  solver_ = continuousSolver_; 

4341  // move solution to continuous copy 

4342  solver_>setColSolution(solution); 

4343  // Put current solution in safe place 

4344  // Point to current solution 

4345  const double * save = testSolution_; 

4346  // Safe as will be const inside infeasibility() 

4347  testSolution_ = solver_>getColSolution(); 

4348  //memcpy(currentSolution_,solver_>getColSolution(), 

4349  // numberColumns*sizeof(double)); 

4350  //solver_>messageHandler()>setLogLevel(4); 

4351  

4352  // save original bounds 

4353  double * saveUpper = new double[numberColumns]; 

4354  double * saveLower = new double[numberColumns]; 

4355  memcpy(saveUpper,getColUpper(),numberColumns*sizeof(double)); 

4356  memcpy(saveLower,getColLower(),numberColumns*sizeof(double)); 

4357  

4358  /* 

4359  Run through the objects and use feasibleRegion() to set variable bounds 

4360  so as to fix the variables specified in the objects at their value in this 

4361  solution. Since the object list contains (at least) one object for every 

4362  integer variable, this has the effect of fixing all integer variables. 

4363  */ 

4364  int i; 

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

4366  object_[i]>feasibleRegion(); 

4367  // We can switch off check 

4368  if ((specialOptions_&4)==0) { 

4369  if ((specialOptions_&2)==0) { 

4370  /* 

4371  Remove any existing warm start information to be sure there is no 

4372  residual influence on initialSolve(). 

4373  */ 

4374  CoinWarmStartBasis *slack = 

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

4376  solver_>setWarmStart(slack); 

4377  delete slack ; 

4378  } 

4379  // Give a hint not to do scaling 

4380  //bool saveTakeHint; 

4381  //OsiHintStrength saveStrength; 

4382  //bool gotHint = (solver_>getHintParam(OsiDoScale,saveTakeHint,saveStrength)); 

4383  //assert (gotHint); 

4384  //solver_>setHintParam(OsiDoScale,false,OsiHintTry); 

4385  solver_>initialSolve(); 

4386  //solver_>setHintParam(OsiDoScale,saveTakeHint,saveStrength); 

4387  if (!solver_>isProvenOptimal()) 

4388  { printf("checkSolution infeas! Retrying wihout scaling.\n"); 

4389  bool saveTakeHint; 

4390  OsiHintStrength saveStrength; 

4391  bool savePrintHint; 

4392  solver_>writeMps("infeas"); 

4393  bool gotHint = (solver_>getHintParam(OsiDoReducePrint,savePrintHint,saveStrength)); 

4394  gotHint = (solver_>getHintParam(OsiDoScale,saveTakeHint,saveStrength)); 

4395  solver_>setHintParam(OsiDoScale,false,OsiHintTry); 

4396  solver_>setHintParam(OsiDoReducePrint,false,OsiHintTry) ; 

4397  solver_>initialSolve(); 

4398  solver_>setHintParam(OsiDoScale,saveTakeHint,saveStrength); 

4399  solver_>setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ; 

4400  } 

4401  //assert(solver_>isProvenOptimal()); 

4402  } 

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

4404  

4405  /* 

4406  Check that the solution still beats the objective cutoff. 

4407  

4408  If it passes, make a copy of the primal variable values and do some 

4409  cleanup and checks: 

4410  + Values of all variables are are within original bounds and values of 

4411  all integer variables are within tolerance of integral. 

4412  + There are no constraint violations. 

4413  There really should be no need for the check against original bounds. 

4414  Perhaps an opportunity for a sanity check? 

4415  */ 

4416  if ((solver_>isProvenOptimal()(specialOptions_&4)!=0) && objectiveValue <= cutoff) 

4417  { 

4418  double * solution = new double[numberColumns]; 

4419  memcpy(solution ,solver_>getColSolution(),numberColumns*sizeof(double)) ; 

4420  

4421  const double * rowLower = solver_>getRowLower() ; 

4422  const double * rowUpper = solver_>getRowUpper() ; 

4423  int numberRows = solver_>getNumRows() ; 

4424  double *rowActivity = new double[numberRows] ; 

4425  memset(rowActivity,0,numberRows*sizeof(double)) ; 

4426  

4427  double integerTolerance = getIntegerTolerance() ; 

4428  

4429  int iColumn; 

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

4431  { double value = solution[iColumn] ; 

4432  value = CoinMax(value, saveLower[iColumn]) ; 

4433  value = CoinMin(value, saveUpper[iColumn]) ; 

4434  if (solver_>isInteger(iColumn)) 

4435  assert(fabs(valuesolution[iColumn]) <= integerTolerance) ; 

4436  solution[iColumn] = value ; } 

4437  

4438  solver_>getMatrixByCol()>times(solution,rowActivity) ; 

4439  delete [] solution; 

4440  double primalTolerance ; 

4441  solver_>getDblParam(OsiPrimalTolerance,primalTolerance) ; 

4442  double largestInfeasibility =0.0; 

4443  for (i=0 ; i < numberRows ; i++) { 

4444  largestInfeasibility = CoinMax(largestInfeasibility, 

4445  rowLower[i]rowActivity[i]); 

4446  largestInfeasibility = CoinMax(largestInfeasibility, 

4447  rowActivity[i]rowUpper[i]); 

4448  } 

4449  if (largestInfeasibility>100.0*primalTolerance) 

4450  handler_>message(CBC_NOTFEAS3, messages_) 

4451  << largestInfeasibility << CoinMessageEol ; 

4452  

4453  delete [] rowActivity ; } 

4454  else 

4455  { objectiveValue=1.0e50 ; } 

4456  /* 

4457  Regardless of what we think of the solution, we may need to restore the 

4458  original bounds of the continuous solver. Unfortunately, const'ness 

4459  prevents us from simply reversing the memcpy used to make these snapshots. 

4460  */ 

4461  if (!fixVariables) 

4462  { for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) 

4463  { solver_>setColLower(iColumn,saveLower[iColumn]) ; 

4464  solver_>setColUpper(iColumn,saveUpper[iColumn]) ; } } 

4465  delete [] saveLower; 

4466  delete [] saveUpper; 

4467  

4468  /* 

4469  Restore the usual solver. 

4470  */ 

4471  solver_=saveSolver; 

4472  testSolution_ = save; 

4473  return objectiveValue; 

4474  } 

4475  

4476  /* 

4477  Call this routine from anywhere when a solution is found. The solution 

4478  vector is assumed to contain one value for each structural variable. 

4479  

4480  The first action is to run checkSolution() to confirm the objective and 

4481  feasibility. If this check causes the solution to be rejected, we're done. 

4482  If fixVariables = true, the variable bounds held by the continuous solver 

4483  will be left fixed to the values in the solution; otherwise they are 

4484  restored to the original values. 

4485  

4486  If the solution is accepted, install it as the best solution. 

4487  

4488  The routine also contains a hook to run any cut generators that are flagged 

4489  to run when a new solution is discovered. There's a potential hazard because 

4490  the cut generators see the continuous solver >after< possible restoration of 

4491  original bounds (which may well invalidate the solution). 

4492  */ 

4493  

4494  void 

4495  CbcModel::setBestSolution (CBC_Message how, 

4496  double & objectiveValue, const double *solution, 

4497  bool fixVariables) 

4498  

4499  { double cutoff = getCutoff() ; 

4500  

4501  /* 

4502  Double check the solution to catch pretenders. 

4503  */ 

4504  objectiveValue = checkSolution(cutoff,solution,fixVariables); 

4505  if (objectiveValue > cutoff) 

4506  { if (objectiveValue>1.0e30) 

4507  handler_>message(CBC_NOTFEAS1, messages_) << CoinMessageEol ; 

4508  else 

4509  handler_>message(CBC_NOTFEAS2, messages_) 

4510  << objectiveValue << cutoff << CoinMessageEol ; } 

4511  /* 

4512  We have a winner. Install it as the new incumbent. 

4513  Bump the objective cutoff value and solution counts. Give the user the 

4514  good news. 

4515  */ 

4516  else 

4517  { bestObjective_ = objectiveValue; 

4518  int numberColumns = solver_>getNumCols(); 

4519  if (!bestSolution_) 

4520  bestSolution_ = new double[numberColumns]; 

4521  CoinCopyN(solution,numberColumns,bestSolution_); 

4522  

4523  cutoff = bestObjective_dblParam_[CbcCutoffIncrement]; 

4524  // This is not correct  that way cutoff can go up if maximization 

4525  //double direction = solver_>getObjSense(); 

4526  //setCutoff(cutoff*direction); 

4527  setCutoff(cutoff); 

4528  

4529  if (how==CBC_ROUNDING) 

4530  numberHeuristicSolutions_++; 

4531  numberSolutions_++; 

4532  if (numberHeuristicSolutions_==numberSolutions_) 

4533  stateOfSearch_ = 1; 

4534  else 

4535  stateOfSearch_ = 2; 

4536  

4537  handler_>message(how,messages_) 

4538  <<bestObjective_<<numberIterations_ 

4539  <<numberNodes_ 

4540  <<CoinMessageEol; 

4541  /* 

4542  Now step through the cut generators and see if any of them are flagged to 

4543  run when a new solution is discovered. Only global cuts are useful. (The 

4544  solution being evaluated may not correspond to the current location in the 

4545  search tree  discovered by heuristic, for example.) 

4546  */ 

4547  OsiCuts theseCuts; 

4548  int i; 

4549  int lastNumberCuts=0; 

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

4551  if (generator_[i]>atSolution()) { 

4552  generator_[i]>generateCuts(theseCuts,true,NULL); 

4553  int numberCuts = theseCuts.sizeRowCuts(); 

4554  for (int j=lastNumberCuts;j<numberCuts;j++) { 

4555  const OsiRowCut * thisCut = theseCuts.rowCutPtr(j); 

4556  if (thisCut>globallyValid()) { 

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

4558  /* As these are global cuts  

4559  a) Always get debugger object 

4560  b) Not fatal error to cutoff optimal (if we have just got optimal) 

4561  */ 

4562  const OsiRowCutDebugger *debugger = solver_>getRowCutDebuggerAlways() ; 

4563  if (debugger) { 

4564  if(debugger>invalidCut(*thisCut)) 

4565  printf("ZZZZ Global cut  cuts off optimal solution!\n"); 

4566  } 

4567  } 

4568  // add to global list 

4569  globalCuts_.insert(*thisCut); 

4570  generator_[i]>incrementNumberCutsInTotal(); 

4571  } 

4572  } 

4573  } 

4574  } 

4575  int numberCuts = theseCuts.sizeColCuts(); 

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

4577  const OsiColCut * thisCut = theseCuts.colCutPtr(i); 

4578  if (thisCut>globallyValid()) { 

4579  // add to global list 

4580  globalCuts_.insert(*thisCut); 

4581  } 

4582  } 

4583  } 

4584  

4585  return ; } 

4586  

4587  

4588  /* Test the current solution for feasibility. 

4589  

4590  Calculate the number of standard integer infeasibilities, then scan the 

4591  remaining objects to see if any of them report infeasibilities. 

4592  

4593  Currently (2003.08) the only object besides SimpleInteger is Clique, hence 

4594  the comments about `odd ones' infeasibilities. 

4595  */ 

4596  bool 

4597  CbcModel::feasibleSolution(int & numberIntegerInfeasibilities, 

4598  int & numberObjectInfeasibilities) const 

4599  { 

4600  int numberUnsatisfied=0; 

4601  double sumUnsatisfied=0.0; 

4602  int preferredWay; 

4603  int j; 

4604  // Point to current solution 

4605  const double * save = testSolution_; 

4606  // Safe as will be const inside infeasibility() 

4607  testSolution_ = solver_>getColSolution(); 

4608  // Put current solution in safe place 

4609  //memcpy(currentSolution_,solver_>getColSolution(), 

4610  // solver_>getNumCols()*sizeof(double)); 

4611  for (j=0;j<numberIntegers_;j++) { 

4612  const CbcObject * object = object_[j]; 

4613  double infeasibility = object>infeasibility(preferredWay); 

4614  if (infeasibility) { 

4615  assert (infeasibility>0); 

4616  numberUnsatisfied++; 

4617  sumUnsatisfied += infeasibility; 

4618  } 

4619  } 

4620  numberIntegerInfeasibilities = numberUnsatisfied; 

4621  for (;j<numberObjects_;j++) { 

4622  const CbcObject * object = object_[j]; 

4623  double infeasibility = object>infeasibility(preferredWay); 

4624  if (infeasibility) { 

4625  assert (infeasibility>0); 

4626  numberUnsatisfied++; 

4627  sumUnsatisfied += infeasibility; 

4628  } 

4629  } 

4630  // and restore 

4631  testSolution_ = save; 

4632  numberObjectInfeasibilities = numberUnsatisfiednumberIntegerInfeasibilities; 

4633  return (!numberUnsatisfied); 

4634  } 

4635  

4636  /* For all vubs see if we can tighten bounds by solving Lp's 

4637  type  0 just vubs 

4638  1 all (could be very slow) 

4639  1 just vubs where variable away from bound 

4640  Returns false if not feasible 

4641  */ 

4642  bool 

4643  CbcModel::tightenVubs(int type, bool allowMultipleBinary, double useCutoff) 

4644  { 

4645  

4646  CoinPackedMatrix matrixByRow(*solver_>getMatrixByRow()); 

4647  int numberRows = solver_>getNumRows(); 

4648  int numberColumns = solver_>getNumCols(); 

4649  

4650  int iRow,iColumn; 

4651  

4652  // Row copy 

4653  //const double * elementByRow = matrixByRow.getElements(); 

4654  const int * column = matrixByRow.getIndices(); 

4655  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts(); 

4656  const int * rowLength = matrixByRow.getVectorLengths(); 

4657  

4658  const double * colUpper = solver_>getColUpper(); 

4659  const double * colLower = solver_>getColLower(); 

4660  //const double * rowUpper = solver_>getRowUpper(); 

4661  //const double * rowLower = solver_>getRowLower(); 

4662  

4663  const double * objective = solver_>getObjCoefficients(); 

4664  //double direction = solver_>getObjSense(); 

4665  const double * colsol = solver_>getColSolution(); 

4666  

4667  int numberVub=0; 

4668  int * continuous = new int[numberColumns]; 

4669  if (type >= 0) { 

4670  double * sort = new double[numberColumns]; 

4671  for (iRow=0;iRow<numberRows;iRow++) { 

4672  int j; 

4673  int numberBinary=0; 

4674  int numberUnsatisfiedBinary=0; 

4675  int numberContinuous=0; 

4676  int iCont=1; 

4677  double weight=1.0e30; 

4678  for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) { 

4679  int iColumn = column[j]; 

4680  if (colUpper[iColumn]colLower[iColumn]>1.0e8) { 

4681  if (solver_>isFreeBinary(iColumn)) { 

4682  numberBinary++; 

4683  /* For sort I make naive assumption: 

4684  x  a * delta <=0 or 

4685  x + a * delta >= 0 

4686  */ 

4687  if (colsol[iColumn]>colLower[iColumn]+1.0e6&& 

4688  colsol[iColumn]<colUpper[iColumn]1.0e6) { 

4689  numberUnsatisfiedBinary++; 

4690  weight = CoinMin(weight,fabs(objective[iColumn])); 

4691  } 

4692  } else { 

4693  numberContinuous++; 

4694  iCont=iColumn; 

4695  } 

4696  } 

4697  } 

4698  if (numberContinuous==1&&numberBinary) { 

4699  if (numberBinary==1allowMultipleBinary) { 

4700  // treat as vub 

4701  if (!numberUnsatisfiedBinary) 

4702  weight=1.0; // at end 

4703  sort[numberVub]=weight; 

4704  continuous[numberVub++] = iCont; 

4705  } 

4706  } 

4707  } 

4708  if (type>0) { 

4709  // take so many 

4710  CoinSort_2(sort,sort+numberVub,continuous); 

4711  numberVub = CoinMin(numberVub,type); 

4712  } 

4713  delete [] sort; 

4714  } else { 

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

4716  continuous[iColumn]=iColumn; 

4717  numberVub=numberColumns; 

4718  } 

4719  bool feasible = tightenVubs(numberVub,continuous,useCutoff); 

4720  delete [] continuous; 

4721  

4722  return feasible; 

4723  } 

4724  // This version is just handed a list of variables 

4725  bool 

4726  CbcModel::tightenVubs(int numberSolves, const int * which, 

4727  double useCutoff) 

4728  { 

4729  

4730  int numberColumns = solver_>getNumCols(); 

4731  

4732  int iColumn; 

4733  

4734  OsiSolverInterface * solver = solver_; 

4735  double saveCutoff = getCutoff() ; 

4736  

4737  double * objective = new double[numberColumns]; 

4738  memcpy(objective,solver_>getObjCoefficients(),numberColumns*sizeof(double)); 

4739  double direction = solver_>getObjSense(); 

4740  

4741  // add in objective if there is a cutoff 

4742  if (useCutoff<1.0e30) { 

4743  // get new version of model 

4744  solver = solver_>clone(); 

4745  CoinPackedVector newRow; 

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

4747  solver>setObjCoeff(iColumn,0.0); // zero out in new model 

4748  if (objective[iColumn]) 

4749  newRow.insert(iColumn,direction * objective[iColumn]); 

4750  

4751  } 

4752  solver>addRow(newRow,COIN_DBL_MAX,useCutoff); 

4753  // signal no objective 

4754  delete [] objective; 

4755  objective=NULL; 

4756  } 

4757  setCutoff(COIN_DBL_MAX); 

4758  

4759  

4760  bool * vub = new bool [numberColumns]; 

4761  int iVub; 

4762  

4763  // mark vub columns 

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

4765  vub[iColumn]=false; 

4766  for (iVub=0;iVub<numberSolves;iVub++) 

4767  vub[which[iVub]]=true; 

4768  OsiCuts cuts; 

4769  // First tighten bounds anyway if CglProbing there 

4770  CglProbing* generator = NULL; 

4771  int iGen; 

4772  for (iGen=0;iGen<numberCutGenerators_;iGen++) { 

4773  generator = dynamic_cast<CglProbing*>(generator_[iGen]>generator()); 

4774  if (generator) 

4775  break; 

4776  } 

4777  int numberFixed=0; 

4778  int numberTightened=0; 

4779  int numberFixedByProbing=0; 

4780  int numberTightenedByProbing=0; 

4781  int printFrequency = (numberSolves+19)/20; // up to 20 messages 

4782  int save[4]; 

4783  if (generator) { 

4784  // set to cheaper and then restore at end 

4785  save[0]=generator>getMaxPass(); 

4786  save[1]=generator>getMaxProbe(); 

4787  save[2]=generator>getMaxLook(); 

4788  save[3]=generator>rowCuts(); 

4789  generator>setMaxPass(1); 

4790  generator>setMaxProbe(10); 

4791  generator>setMaxLook(50); 

4792  generator>setRowCuts(0); 

4793  

4794  // Probing  return tight column bounds 

4795  generator>generateCutsAndModify(*solver,cuts); 

4796  const double * tightLower = generator>tightLower(); 

4797  const double * lower = solver>getColLower(); 

4798  const double * tightUpper = generator>tightUpper(); 

4799  const double * upper = solver>getColUpper(); 

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

4801  double newUpper = tightUpper[iColumn]; 

4802  double newLower = tightLower[iColumn]; 

4803  if (newUpper<upper[iColumn]1.0e8*(fabs(upper[iColumn])+1) 

4804  newLower>lower[iColumn]+1.0e8*(fabs(lower[iColumn])+1)) { 

4805  if (newUpper<newLower) { 

4806  fprintf(stderr,"Problem is infeasible\n"); 

4807  return false; 

4808  } 

4809  if (newUpper==newLower) { 

4810  numberFixed++; 

4811  numberFixedByProbing++; 

4812  solver>setColLower(iColumn,newLower); 

4813  solver>setColUpper(iColumn,newUpper); 

4814  printf("Column %d, new bounds %g %g\n",iColumn, 

4815  newLower,newUpper); 

4816  } else if (vub[iColumn]) { 

4817  numberTightened++; 

4818  numberTightenedByProbing++; 

4819  if (!solver>isInteger(iColumn)) { 

4820  // relax 

4821  newLower=CoinMax(lower[iColumn], 

4822  newLower 

4823  1.0e5*(fabs(lower[iColumn])+1)); 

4824  newUpper=CoinMin(upper[iColumn], 

4825  newUpper 

4826  +1.0e5*(fabs(upper[iColumn])+1)); 

4827  } 

4828  solver>setColLower(iColumn,newLower); 

4829  solver>setColUpper(iColumn,newUpper); 

4830  } 

4831  } 

4832  } 

4833  } 

4834  CoinWarmStart * ws = solver>getWarmStart(); 

4835  double * solution = new double [numberColumns]; 

4836  memcpy(solution,solver>getColSolution(),numberColumns*sizeof(double)); 

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

4838  solver>setObjCoeff(iColumn,0.0); 

4839  //solver>messageHandler()>setLogLevel(2); 

4840  for (iVub=0;iVub<numberSolves;iVub++) { 

4841  iColumn = which[iVub]; 

4842  int iTry; 

4843  for (iTry=0;iTry<2;iTry++) { 

4844  double saveUpper = solver>getColUpper()[iColumn]; 

4845  double saveLower = solver>getColLower()[iColumn]; 

4846  double value; 

4847  if (iTry==1) { 

4848  // try all way up 

4849  solver>setObjCoeff(iColumn,1.0); 

4850  } else { 

4851  // try all way down 

4852  solver>setObjCoeff(iColumn,1.0); 

4853  } 

4854  solver>initialSolve(); 

4855  value = solver>getColSolution()[iColumn]; 

4856  bool change=false; 

4857  if (iTry==1) { 

4858  if (value<saveUpper1.0e4) { 

4859  if (solver>isInteger(iColumn)) { 

4860  value = floor(value+0.00001); 

4861  } else { 

4862  // relax a bit 

4863  value=CoinMin(saveUpper,value+1.0e5*(fabs(saveUpper)+1)); 

4864  } 

4865  if (valuesaveLower<1.0e7) 

4866  value = saveLower; // make sure exactly same 

4867  solver>setColUpper(iColumn,value); 

4868  saveUpper=value; 

4869  change=true; 

4870  } 

4871  } else { 

4872  if (value>saveLower+1.0e4) { 

4873  if (solver>isInteger(iColumn)) { 

4874  value = ceil(value0.00001); 

4875  } else { 

4876  // relax a bit 

4877  value=CoinMax(saveLower,value1.0e5*(fabs(saveLower)+1)); 

4878  } 

4879  if (saveUppervalue<1.0e7) 

4880  value = saveUpper; // make sure exactly same 

4881  solver>setColLower(iColumn,value); 

4882  saveLower=value; 

4883  change=true; 

4884  } 

4885  } 

4886  solver>setObjCoeff(iColumn,0.0); 

4887  if (change) { 

4888  if (saveUpper==saveLower) 

4889  numberFixed++; 

4890  else 

4891  numberTightened++; 

4892  int saveFixed=numberFixed; 

4893  

4894  int jColumn; 

4895  if (generator) { 

4896  // Probing  return tight column bounds 

4897  cuts = OsiCuts(); 

4898  generator>generateCutsAndModify(*solver,cuts); 

4899  const double * tightLower = generator>tightLower(); 

4900  const double * lower = solver>getColLower(); 

4901  const double * tightUpper = generator>tightUpper(); 

4902  const double * upper = solver>getColUpper(); 

4903  for (jColumn=0;jColumn<numberColumns;jColumn++) { 

4904  double newUpper = tightUpper[jColumn]; 

4905  double newLower = tightLower[jColumn]; 

4906  if (newUpper<upper[jColumn]1.0e8*(fabs(upper[jColumn])+1) 

4907  newLower>lower[jColumn]+1.0e8*(fabs(lower[jColumn])+1)) { 

4908  if (newUpper<newLower) { 

4909  fprintf(stderr,"Problem is infeasible\n"); 

4910  return false; 

4911  } 

4912  if (newUpper==newLower) { 

4913  numberFixed++; 

4914  numberFixedByProbing++; 

4915  solver>setColLower(jColumn,newLower); 

4916  solver>setColUpper(jColumn,newUpper); 

4917  } else if (vub[jColumn]) { 

4918  numberTightened++; 

4919  numberTightenedByProbing++; 

4920  if (!solver>isInteger(jColumn)) { 

4921  // relax 

4922  newLower=CoinMax(lower[jColumn], 

4923  newLower 

4924  1.0e5*(fabs(lower[jColumn])+1)); 

4925  newUpper=CoinMin(upper[jColumn], 

4926  newUpper 

4927  +1.0e5*(fabs(upper[jColumn])+1)); 

4928  } 

4929  solver>setColLower(jColumn,newLower); 

4930  solver>setColUpper(jColumn,newUpper); 

4931  } 

4932  } 

4933  } 

4934  } 

4935  if (numberFixed>saveFixed) { 

4936  // original solution may not be feasible 

4937  // go back to true costs to solve if exists 

4938  if (objective) { 

4939  for (jColumn=0;jColumn<numberColumns;jColumn++) 

4940  solver>setObjCoeff(jColumn,objective[jColumn]); 

4941  } 

4942  solver>setColSolution(solution); 

4943  solver>setWarmStart(ws); 

4944  solver>resolve(); 

4945  if (!solver>isProvenOptimal()) { 

4946  fprintf(stderr,"Problem is infeasible\n"); 

4947  return false; 

4948  } 

4949  delete ws; 

4950  ws = solver>getWarmStart(); 

4951  memcpy(solution,solver>getColSolution(), 

4952  numberColumns*sizeof(double)); 

4953  for (jColumn=0;jColumn<numberColumns;jColumn++) 

4954  solver>setObjCoeff(jColumn,0.0); 

4955  } 

4956  } 

4957  solver>setColSolution(solution); 

4958  solver>setWarmStart(ws); 

4959  } 

4960  if (iVub%printFrequency==0) 

4961  handler_>message(CBC_VUB_PASS,messages_) 

4962  <<iVub+1<<numberFixed<<numberTightened 

4963  <<CoinMessageEol; 

4964  } 

4965  handler_>message(CBC_VUB_END,messages_) 

4966  <<numberFixed<<numberTightened 

4967  <<CoinMessageEol; 

4968  delete ws; 

4969  delete [] solution; 

4970  // go back to true costs to solve if exists 

4971  if (objective) { 

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

4973  solver_>setObjCoeff(iColumn,objective[iColumn]); 

4974  delete [] objective; 

4975  } 

4976  delete [] vub; 

4977  if (generator) { 

4978  /*printf("Probing fixed %d and tightened %d\n", 

4979  numberFixedByProbing, 

4980  numberTightenedByProbing);*/ 

4981  if (generator_[iGen]>howOften()==1&& 

4982  (numberFixedByProbing+numberTightenedByProbing)*5> 

4983  (numberFixed+numberTightened)) 

4984  generator_[iGen]>setHowOften(1000000+1); 

4985  generator>setMaxPass(save[0]); 

4986  generator>setMaxProbe(save[1]); 

4987  generator>setMaxLook(save[2]); 

4988  generator>setRowCuts(save[3]); 

4989  } 

4990  

4991  if (solver!=solver_) { 

4992  // move bounds across 

4993  const double * lower = solver>getColLower(); 

4994  const double * upper = solver>getColUpper(); 

4995  const double * lowerOrig = solver_>getColLower(); 

4996  const double * upperOrig = solver_>getColUpper(); 

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

4998  solver_>setColLower(iColumn,CoinMax(lower[iColumn],lowerOrig[iColumn])); 

4999  solver_>setColUpper(iColumn,CoinMin(upper[iColumn],upperOrig[iColumn])); 

5000  } 

5001  delete solver; 

5002  } 

5003  setCutoff(saveCutoff); 

5004  return true; 

5005  } 

5006  /* 

5007  Do Integer Presolve. Returns new model. 

5008  I have to work out cleanest way of getting solution to 

5009  original problem at end. So this is very preliminary. 

5010  */ 

5011  CbcModel * 

5012  CbcModel::integerPresolve(bool weak) 

5013  { 

5014  status_ = 0; 

5015  // solve LP 

5016  //solver_>writeMps("bad"); 

5017  bool feasible = resolve(); 

5018  

5019  CbcModel * newModel = NULL; 

5020  if (feasible) { 

5021  

5022  // get a new model 

5023  newModel = new CbcModel(*this); 

5024  newModel>messageHandler()>setLogLevel(messageHandler()>logLevel()); 

5025  

5026  feasible = newModel>integerPresolveThisModel(solver_,weak); 

5027  } 

5028  if (!feasible) { 

5029  handler_>message(CBC_INFEAS,messages_) 

5030  <<CoinMessageEol; 

5031  status_ = 0; 

5032  secondaryStatus_ = 1; 

5033  delete newModel; 

5034  return NULL; 

5035  } else { 

5036  newModel>synchronizeModel(); // make sure everything that needs solver has it 

5037  return newModel; 

5038  } 

5039  } 

5040  /* 

5041  Do Integer Presolve  destroying current model 

5042  */ 

5043  bool 

5044  CbcModel::integerPresolveThisModel(OsiSolverInterface * originalSolver, 

5045  bool weak) 

5046  { 

5047  status_ = 0; 

5048  // solve LP 

5049  bool feasible = resolve(); 

5050  

5051  bestObjective_=1.0e50; 

5052  numberSolutions_=0; 

5053  numberHeuristicSolutions_=0; 

5054  double cutoff = getCutoff() ; 

5055  double direction = solver_>getObjSense(); 

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

5057  messageHandler()>message(CBC_CUTOFF_WARNING1, 

5058  messages()) 

5059  << cutoff << cutoff << CoinMessageEol ; 

5060  if (cutoff > bestObjective_) 

5061  cutoff = bestObjective_ ; 

5062  setCutoff(cutoff) ; 

5063  int iColumn; 

5064  int numberColumns = getNumCols(); 

5065  int originalNumberColumns = numberColumns; 

5066  currentPassNumber_=0; 

5067  synchronizeModel(); // make sure everything that needs solver has it 

5068  // just point to solver_ 

5069  delete continuousSolver_; 

5070  continuousSolver_ = solver_; 

5071  // get a copy of original so we can fix bounds 

5072  OsiSolverInterface * cleanModel = originalSolver>clone(); 

5073  #ifdef CBC_DEBUG 

5074  std::string problemName; 

5075  cleanModel>getStrParam(OsiProbName,problemName); 

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

5077  cleanModel>activateRowCutDebugger(problemName.c_str()); 

5078  const OsiRowCutDebugger * debugger = cleanModel>getRowCutDebugger(); 

5079  #endif 

5080  

5081  // array which points from original columns to presolved 

5082  int * original = new int[numberColumns]; 

5083  // arrays giving bounds  only ones found by probing 

5084  // rest will be found by presolve 

5085  double * originalLower = new double[numberColumns]; 

5086  double * originalUpper = new double[numberColumns]; 

5087  { 

5088  const double * lower = getColLower(); 

5089  const double * upper = getColUpper(); 

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

5091  original[iColumn]=iColumn; 

5092  originalLower[iColumn] = lower[iColumn]; 

5093  originalUpper[iColumn] = upper[iColumn]; 

5094  } 

5095  } 

5096  findIntegers(true); 

5097  // save original integers 

5098  int * originalIntegers = new int[numberIntegers_]; 

5099  int originalNumberIntegers = numberIntegers_; 

5100  memcpy(originalIntegers,integerVariable_,numberIntegers_*sizeof(int)); 

5101  

5102  int todo=20; 

5103  if (weak) 

5104  todo=1; 

5105  while (currentPassNumber_<todo) { 

5106  

5107  currentPassNumber_++; 

5108  numberSolutions_=0; 

5109  // this will be set false to break out of loop with presolved problem 

5110  bool doIntegerPresolve=(currentPassNumber_!=20); 

5111  

5112  // Current number of free integer variables 

5113  // Get increment in solutions 

5114  { 

5115  const double * objective = cleanModel>getObjCoefficients(); 

5116  const double * lower = cleanModel>getColLower(); 

5117  const double * upper = cleanModel>getColUpper(); 

5118  double maximumCost=0.0; 

5119  bool possibleMultiple=true; 

5120  int numberChanged=0; 

5121  for (iColumn=0;iColumn<originalNumberColumns;iColumn++) { 

5122  if (originalUpper[iColumn]>originalLower[iColumn]) { 

5123  if( cleanModel>isInteger(iColumn)) { 

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

5125  } else if (objective[iColumn]) { 

5126  possibleMultiple=false; 

5127  } 

5128  } 

5129  if (originalUpper[iColumn]<upper[iColumn]) { 

5130  #ifdef CBC_DEBUG 

5131  printf("Changing upper bound on %d from %g to %g\n", 

5132  iColumn,upper[iColumn],originalUpper[iColumn]); 

5133  #endif 

5134  cleanModel>setColUpper(iColumn,originalUpper[iColumn]); 

5135  numberChanged++; 

5136  } 

5137  if (originalLower[iColumn]>lower[iColumn]) { 

5138  #ifdef CBC_DEBUG 

5139  printf("Changing lower bound on %d from %g to %g\n", 

5140  iColumn,lower[iColumn],originalLower[iColumn]); 

5141  #endif 

5142  cleanModel>setColLower(iColumn,originalLower[iColumn]); 

5143  numberChanged++; 

5144  } 

5145  } 

5146  // if first pass  always try 

5147  if (currentPassNumber_==1) 

5148  numberChanged += 1; 

5149  if (possibleMultiple&&maximumCost) { 

5150  int increment=0; 

5151  double multiplier = 2520.0; 

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

5153  multiplier *= 10.0; 

5154  for (int j =0;j<originalNumberIntegers;j++) { 

5155  iColumn = originalIntegers[j]; 

5156  if (originalUpper[iColumn]>originalLower[iColumn]) { 

5157  if(objective[iColumn]) { 

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

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

5160  if (fabs(valuefloor(value+0.5))>1.0e8) { 

5161  increment=0; 

5162  break; // no good 

5163  } else if (!increment) { 

5164  // first 

5165  increment=nearest; 

5166  } else { 

5167  increment = gcd(increment,nearest); 

5168  } 

5169  } 

5170  } 

5171  } 

5172  if (increment) { 

5173  double value = increment; 

5174  value /= multiplier; 

5175  if (value*0.999>dblParam_[CbcCutoffIncrement]) { 

5176  messageHandler()>message(CBC_INTEGERINCREMENT,messages()) 

5177  <<value 

5178  <<CoinMessageEol; 

5179  dblParam_[CbcCutoffIncrement]=value*0.999; 

5180  } 

5181  } 

5182  } 

5183  if (!numberChanged) { 

5184  doIntegerPresolve=false; // not doing any better 

5185  } 

5186  } 

5187  #ifdef CBC_DEBUG 

5188  if (debugger) 

5189  assert(debugger>onOptimalPath(*cleanModel)); 

5190  #endif 

5191  #ifdef COIN_USE_CLP 

5192  // do presolve  for now just clp but easy to get osi interface 

5193  OsiClpSolverInterface * clpSolver 

5194  = dynamic_cast<OsiClpSolverInterface *> (cleanModel); 

5195  if (clpSolver) { 

5196  ClpSimplex * clp = clpSolver>getModelPtr(); 

5197  clp>messageHandler()>setLogLevel(cleanModel>messageHandler()>logLevel()); 

5198  ClpPresolve pinfo; 

5199  //printf("integerPresolve  temp switch off doubletons\n"); 

5200  //pinfo.setPresolveActions(4); 

5201  ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e8); 

5202  if (!model2) { 

5203  // presolve found to be infeasible 

5204  feasible=false; 

5205  } else { 

5206  // update original array 

5207  const int * originalColumns = pinfo.originalColumns(); 

5208  // just slot in new solver 

5209  OsiClpSolverInterface * temp = new OsiClpSolverInterface(model2,true); 

5210  numberColumns = temp>getNumCols(); 

5211  for (iColumn=0;iColumn<originalNumberColumns;iColumn++) 

5212  original[iColumn]=1; 

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

5214  original[originalColumns[iColumn]]=iColumn; 

5215  // copy parameters 

5216  temp>copyParameters(*solver_); 

5217  // and specialized ones 

5218  temp>setSpecialOptions(clpSolver>specialOptions()); 

5219  delete solver_; 

5220  solver_ = temp; 

5221  setCutoff(cutoff); 

5222  deleteObjects(); 

5223  if (!numberObjects_) { 

5224  // Nothing left 

5225  doIntegerPresolve=false; 

5226  weak=true; 

5227  break; 

5228  } 

5229  synchronizeModel(); // make sure everything that needs solver has it 

5230  // just point to solver_ 

5231  continuousSolver_ = solver_; 

5232  feasible=resolve(); 

5233  if (!feasible!doIntegerPresolveweak) break; 

5234  // see if we can get solution by heuristics 

5235  int found=1; 

5236  int iHeuristic; 

5237  double * newSolution = new double [numberColumns]; 

5238  double heuristicValue=getCutoff(); 

5239  for (iHeuristic=0;iHeuristic<numberHeuristics_;iHeuristic++) { 

5240  double saveValue=heuristicValue; 

5241  int ifSol = heuristic_[iHeuristic]>solution(heuristicValue, 

5242  newSolution); 

5243  if (ifSol>0) { 

5244  // better solution found 

5245  found=iHeuristic; 

5246  incrementUsed(newSolution); 

5247  } else if (ifSol<0) { 

5248  heuristicValue = saveValue; 

5249  } 

5250  } 

5251  if (found >= 0) { 

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

5253  int numberColumns = getNumCols() ; 

5254  if (!currentSolution_) 

5255  currentSolution_ = new double[numberColumns] ; 

5256  testSolution_=currentSolution_; 

5257  // better solution save 

5258  setBestSolution(CBC_ROUNDING,heuristicValue, 

5259  newSolution); 

5260  // update cutoff 

5261  cutoff = getCutoff(); 

5262  } 

5263  delete [] newSolution; 

5264  // Space for type of cuts 

5265  int maximumWhich=1000; 

5266  int * whichGenerator = new int[maximumWhich]; 

5267  // save number of rows 

5268  numberRowsAtContinuous_ = getNumRows(); 

5269  maximumNumberCuts_=0; 

5270  currentNumberCuts_=0; 

5271  delete [] addedCuts_; 

5272  addedCuts_ = NULL; 

5273  

5274  // maximum depth for tree walkback 

5275  maximumDepth_=10; 

5276  delete [] walkback_; 

5277  walkback_ = new CbcNodeInfo * [maximumDepth_]; 

5278  

5279  OsiCuts cuts; 

5280  int numberOldActiveCuts=0; 

5281  int numberNewCuts = 0; 

5282  feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_, 

5283  NULL,numberOldActiveCuts,numberNewCuts, 

5284  maximumWhich, whichGenerator); 

5285  currentNumberCuts_=numberNewCuts; 

5286  delete [] whichGenerator; 

5287  delete [] walkback_; 

5288  walkback_ = NULL; 

5289  delete [] addedCuts_; 

5290  addedCuts_=NULL; 

5291  if (feasible) { 

5292  // fix anything in original which integer presolve fixed 

5293  // for now just integers 

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

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

5296  int i; 

5297  for (i=0;i<originalNumberIntegers;i++) { 

5298  iColumn = originalIntegers[i]; 

5299  int jColumn = original[iColumn]; 

5300  if (jColumn >= 0) { 

5301  if (upper[jColumn]<originalUpper[iColumn]) 

5302  originalUpper[iColumn] = upper[jColumn]; 

5303  if (lower[jColumn]>originalLower[iColumn]) 

5304  originalLower[iColumn] = lower[jColumn]; 

5305  } 

5306  } 

5307  } 

5308  } 

5309  } 

5310  #endif 

5311  if (!feasible!doIntegerPresolve) { 

5312  break; 

5313  } 

5314  } 

5315  //solver_>writeMps("xx"); 

5316  delete cleanModel; 

5317  delete [] originalIntegers; 

5318  numberColumns = getNumCols(); 

5319  delete [] originalColumns_; 

5320  originalColumns_ = new int[numberColumns]; 

5321  numberColumns=0; 

5322  for (iColumn=0;iColumn<originalNumberColumns;iColumn++) { 

5323  int jColumn = original[iColumn]; 

5324  if (jColumn >= 0) 

5325  originalColumns_[numberColumns++]=iColumn; 

5326  } 

5327  delete [] original; 

5328  delete [] originalLower; 

5329  delete [] originalUpper; 

5330  

5331  deleteObjects(); 

5332  synchronizeModel(); // make sure everything that needs solver has it 

5333  continuousSolver_=NULL; 

5334  currentNumberCuts_=0; 

5335  return feasible; 

5336  } 

5337  // Put back information into original model  after integerpresolve 

5338  void 

5339  CbcModel::originalModel(CbcModel * presolvedModel,bool weak) 

5340  { 

5341  solver_>copyParameters(*(presolvedModel>solver_)); 

5342  bestObjective_ = presolvedModel>bestObjective_; 

5343  delete [] bestSolution_; 

5344  findIntegers(true); 

5345  if (presolvedModel>bestSolution_) { 

5346  int numberColumns = getNumCols(); 

5347  int numberOtherColumns = presolvedModel>getNumCols(); 

5348  //bestSolution_ = new double[numberColumns]; 

5349  // set up map 

5350  int * back = new int[numberColumns]; 

5351  int i; 

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

5353  back[i]=1; 

5354  for (i=0;i<numberOtherColumns;i++) 

5355  back[presolvedModel>originalColumns_[i]]=i; 

5356  int iColumn; 

5357  // set ones in presolved model to values 

5358  double * otherSolution = presolvedModel>bestSolution_; 

5359  //const double * lower = getColLower(); 

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

5361  iColumn = integerVariable_[i]; 

5362  int jColumn = back[iColumn]; 

5363  //bestSolution_[iColumn]=lower[iColumn]; 

5364  if (jColumn >= 0) { 

5365  double value=floor(otherSolution[jColumn]+0.5); 

5366  solver_>setColLower(iColumn,value); 

5367  solver_>setColUpper(iColumn,value); 

5368  //bestSolution_[iColumn]=value; 

5369  } 

5370  } 

5371  delete [] back; 

5372  #if 0 

5373  // ** looks as if presolve needs more intelligence 

5374  // do presolve  for now just clp but easy to get osi interface 

5375  OsiClpSolverInterface * clpSolver 

5376  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

5377  assert (clpSolver); 

5378  ClpSimplex * clp = clpSolver>getModelPtr(); 

5379  Presolve pinfo; 

5380  ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e8); 

5381  model2>primal(1); 

5382  pinfo.postsolve(true); 

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

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

5385  iColumn = integerVariable_[i]; 

5386  double value=floor(solution[iColumn]+0.5); 

5387  solver_>setColLower(iColumn,value); 

5388  solver_>setColUpper(iColumn,value); 

5389  } 

5390  #else 

5391  if (!weak) { 

5392  // for now give up 

5393  int save = numberCutGenerators_; 

5394  numberCutGenerators_=0; 

5395  bestObjective_=1.0e100; 

5396  branchAndBound(); 

5397  numberCutGenerators_=save; 

5398  } 

5399  #endif 

5400  if (bestSolution_) { 

5401  // solve problem 

5402  resolve(); 

5403  // should be feasible 

5404  int numberIntegerInfeasibilities; 

5405  int numberObjectInfeasibilities; 

5406  if (!currentSolution_) 

5407  currentSolution_ = new double[numberColumns] ; 

5408  testSolution_ = currentSolution_; 

5409  assert(feasibleSolution(numberIntegerInfeasibilities, 

5410  numberObjectInfeasibilities)); 

5411  } 

5412  } else { 

5413  bestSolution_=NULL; 

5414  } 

5415  numberSolutions_=presolvedModel>numberSolutions_; 

5416  numberHeuristicSolutions_=presolvedModel>numberHeuristicSolutions_; 

5417  numberNodes_ = presolvedModel>numberNodes_; 

5418  numberIterations_ = presolvedModel>numberIterations_; 

5419  status_ = presolvedModel>status_; 

5420  secondaryStatus_ = presolvedModel>secondaryStatus_; 

5421  synchronizeModel(); 

5422  } 

5423  // Pass in Message handler (not deleted at end) 

5424  void 

5425  CbcModel::passInMessageHandler(CoinMessageHandler * handler) 

5426  { 

5427  if (defaultHandler_) { 

5428  delete handler_; 

5429  handler_ = NULL; 

5430  } 

5431  defaultHandler_=false; 

5432  handler_=handler; 

5433  } 

5434  void 

5435  CbcModel::passInTreeHandler(CbcTree & tree) 

5436  { 

5437  delete tree_; 

5438  tree_ = tree.clone(); 

5439  } 

5440  // Make sure region there 

5441  void 

5442  CbcModel::reserveCurrentSolution(const double * solution) 

5443  { 

5444  int numberColumns = getNumCols() ; 

5445  if (!currentSolution_) 

5446  currentSolution_ = new double[numberColumns] ; 

5447  testSolution_=currentSolution_; 

5448  if (solution) 

5449  memcpy(currentSolution_,solution,numberColumns*sizeof(double)); 

5450  } 

5451  /* For passing in an CbcModel to do a sub Tree (with derived tree handlers). 

5452  Passed in model must exist for duration of branch and bound 

5453  */ 

5454  void 

5455  CbcModel::passInSubTreeModel(CbcModel & model) 

5456  { 

5457  subTreeModel_=&model; 

5458  } 

5459  // For retrieving a copy of subtree model with given OsiSolver or NULL 

5460  CbcModel * 

5461  CbcModel::subTreeModel(OsiSolverInterface * solver) const 

5462  { 

5463  const CbcModel * subModel=subTreeModel_; 

5464  if (!subModel) 

5465  subModel=this; 

5466  // Get new copy 

5467  CbcModel * newModel = new CbcModel(*subModel); 

5468  if (solver) 

5469  newModel>assignSolver(solver); 

5470  return newModel; 

5471  } 

5472  //############################################################################# 

5473  // Set/Get Application Data 

5474  // This is a pointer that the application can store into and retrieve 

5475  // from the solverInterface. 

5476  // This field is the application to optionally define and use. 

5477  //############################################################################# 

5478  

5479  void CbcModel::setApplicationData(void * appData) 

5480  { 

5481  appData_ = appData; 

5482  } 

5483  // 

5484  void * CbcModel::getApplicationData() const 

5485  { 

5486  return appData_; 

5487  } 

5488  /* create a submodel from partially fixed problem 

5489  

5490  The method creates a new clean model with given bounds. 

5491  */ 

5492  CbcModel * 

5493  CbcModel::cleanModel(const double * lower, const double * upper) 

5494  { 

5495  OsiSolverInterface * solver = continuousSolver_>clone(); 

5496  

5497  int numberIntegers = numberIntegers_; 

5498  const int * integerVariable = integerVariable_; 

5499  

5500  int i; 

5501  for (i=0;i<numberIntegers;i++) { 

5502  int iColumn=integerVariable[i]; 

5503  const CbcObject * object = object_[i]; 

5504  const CbcSimpleInteger * integerObject = 

5505  dynamic_cast<const CbcSimpleInteger *> (object); 

5506  assert(integerObject); 

5507  // get original bounds 

5508  double originalLower = integerObject>originalLowerBound(); 

5509  double originalUpper = integerObject>originalUpperBound(); 

5510  solver>setColLower(iColumn,CoinMax(lower[iColumn],originalLower)); 

5511  solver>setColUpper(iColumn,CoinMin(upper[iColumn],originalUpper)); 

5512  } 

5513  CbcModel * model = new CbcModel(*solver); 

5514  // off some messages 

5515  if (handler_>logLevel()<=1) { 

5516  model>messagesPointer()>setDetailMessage(3,9); 

5517  model>messagesPointer()>setDetailMessage(3,6); 

5518  model>messagesPointer()>setDetailMessage(3,4); 

5519  model>messagesPointer()>setDetailMessage(3,1); 

5520  model>messagesPointer()>setDetailMessage(3,13); 

5521  model>messagesPointer()>setDetailMessage(3,14); 

5522  model>messagesPointer()>setDetailMessage(3,3007); 

5523  } 

5524  // Cuts 

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

5526  int howOften = generator_[i]>howOftenInSub(); 

5527  if (howOften>100) { 

5528  CbcCutGenerator * generator = virginGenerator_[i]; 

5529  CglCutGenerator * cglGenerator = generator>generator(); 

5530  model>addCutGenerator(cglGenerator,howOften, 

5531  generator>cutGeneratorName(), 

5532  generator>normal(), 

5533  generator>atSolution(), 

5534  generator>whenInfeasible(), 

5535  100, generator>whatDepthInSub(),1); 

5536  } 

5537  } 

5538  double cutoff = getCutoff(); 

5539  model>setCutoff(cutoff); 

5540  return model; 

5541  } 

5542  /* Invoke the branch & cut algorithm on partially fixed problem 

5543  

5544  The method uses a subModel created by cleanModel. The search 

5545  ends when the tree is exhausted or maximum nodes is reached. 

5546  

5547  If better solution found then it is saved. 

5548  

5549  Returns 0 if search completed and solution, 1 if not completed and solution, 

5550  2 if completed and no solution, 3 if not completed and no solution. 

5551  

5552  Normally okay to do subModel immediately followed by subBranchandBound 

5553  (== other form of subBranchAndBound) 

5554  but may need to get at model for advanced features. 

5555  

5556  Deletes model 

5557  

5558  */ 

5559  

5560  int 

5561  CbcModel::subBranchAndBound(CbcModel * model, 

5562  CbcModel * presolvedModel, 

5563  int maximumNodes) 

5564  { 

5565  int i; 

5566  double cutoff=model>getCutoff(); 

5567  CbcModel * model2; 

5568  if (presolvedModel) 

5569  model2=presolvedModel; 

5570  else 

5571  model2=model; 

5572  // Do complete search 

5573  

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

5575  model2>addHeuristic(heuristic_[i]); 

5576  model2>heuristic(i)>resetModel(model2); 

5577  } 

5578  // Definition of node choice 

5579  model2>setNodeComparison(nodeCompare_>clone()); 

5580  //model2>solver()>setHintParam(OsiDoReducePrint,true,OsiHintTry); 

5581  model2>messageHandler()>setLogLevel(CoinMax(0,handler_>logLevel()1)); 

5582  //model2>solver()>messageHandler()>setLogLevel(2); 

5583  model2>setMaximumCutPassesAtRoot(maximumCutPassesAtRoot_); 

5584  model2>setPrintFrequency(50); 

5585  model2>setIntParam(CbcModel::CbcMaxNumNode,maximumNodes); 

5586  model2>branchAndBound(); 

5587  delete model2>nodeComparison(); 

5588  if (model2>getMinimizationObjValue()>cutoff) { 

5589  // no good 

5590  if (model!=model2) 

5591  delete model2; 

5592  delete model; 

5593  return 2; 

5594  } 

5595  if (model!=model2) { 

5596  // get back solution 

5597  model>originalModel(model2,false); 

5598  delete model2; 

5599  } 

5600  int status; 

5601  if (model>getMinimizationObjValue()<cutoff&&model>bestSolution()) { 

5602  double objValue = model>getObjValue(); 

5603  const double * solution = model>bestSolution(); 

5604  setBestSolution(CBC_TREE_SOL,objValue,solution); 

5605  status = 0; 

5606  } else { 

5607  status=2; 

5608  } 

5609  if (model>status()) 

5610  status ++ ; // not finished search 

5611  delete model; 

5612  return status; 

5613  } 

5614  /* Invoke the branch & cut algorithm on partially fixed problem 

5615  

5616  The method creates a new model with given bounds, presolves it 

5617  then proceeds to explore the branch & cut search tree. The search 

5618  ends when the tree is exhausted or maximum nodes is reached. 

5619  Returns 0 if search completed and solution, 1 if not completed and solution, 

5620  2 if completed and no solution, 3 if not completed and no solution. 

5621  */ 

5622  int 

5623  CbcModel::subBranchAndBound(const double * lower, const double * upper, 

5624  int maximumNodes) 

5625  { 

5626  OsiSolverInterface * solver = continuousSolver_>clone(); 

5627  

5628  int numberIntegers = numberIntegers_; 

5629  const int * integerVariable = integerVariable_; 

5630  

5631  int i; 

5632  for (i=0;i<numberIntegers;i++) { 

5633  int iColumn=integerVariable[i]; 

5634  const CbcObject * object = object_[i]; 

5635  const CbcSimpleInteger * integerObject = 

5636  dynamic_cast<const CbcSimpleInteger *> (object); 

5637  assert(integerObject); 

5638  // get original bounds 

5639  double originalLower = integerObject>originalLowerBound(); 

5640  double originalUpper = integerObject>originalUpperBound(); 

5641  solver>setColLower(iColumn,CoinMax(lower[iColumn],originalLower)); 

5642  solver>setColUpper(iColumn,CoinMin(upper[iColumn],originalUpper)); 

5643  } 

5644  CbcModel model(*solver); 

5645  // off some messages 

5646  if (handler_>logLevel()<=1) { 

5647  model.messagesPointer()>setDetailMessage(3,9); 

5648  model.messagesPointer()>setDetailMessage(3,6); 

5649  model.messagesPointer()>setDetailMessage(3,4); 

5650  model.messagesPointer()>setDetailMessage(3,1); 

5651  model.messagesPointer()>setDetailMessage(3,3007); 

5652  } 

5653  double cutoff = getCutoff(); 

5654  model.setCutoff(cutoff); 

5655  // integer presolve 

5656  CbcModel * model2 = model.integerPresolve(false); 

5657  if (!model2!model2>getNumRows()) { 

5658  delete model2; 

5659  delete solver; 

5660  return 2; 

5661  } 

5662  if (handler_>logLevel()>1) 

5663  printf("Reduced model has %d rows and %d columns\n", 

5664  model2>getNumRows(),model2>getNumCols()); 

5665  // Do complete search 

5666  

5667  // Cuts 

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

5669  int howOften = generator_[i]>howOftenInSub(); 

5670  if (howOften>100) { 

5671  CbcCutGenerator * generator = virginGenerator_[i]; 

5672  CglCutGenerator * cglGenerator = generator>generator(); 

5673  model2>addCutGenerator(cglGenerator,howOften, 

5674  generator>cutGeneratorName(), 

5675  generator>normal(), 

5676  generator>atSolution(), 

5677  generator>whenInfeasible(), 

5678  100, generator>whatDepthInSub(),1); 

5679  } 

5680  } 

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

5682  model2>addHeuristic(heuristic_[i]); 

5683  model2>heuristic(i)>resetModel(model2); 

5684  } 

5685  // Definition of node choice 

5686  model2>setNodeComparison(nodeCompare_>clone()); 

5687  //model2>solver()>setHintParam(OsiDoReducePrint,true,OsiHintTry); 

5688  model2>messageHandler()>setLogLevel(CoinMax(0,handler_>logLevel()1)); 

5689  //model2>solver()>messageHandler()>setLogLevel(2); 

5690  model2>setMaximumCutPassesAtRoot(maximumCutPassesAtRoot_); 

5691  model2>setPrintFrequency(50); 

5692  model2>setIntParam(CbcModel::CbcMaxNumNode,maximumNodes); 

5693  model2>branchAndBound(); 

5694  delete model2>nodeComparison(); 

5695  if (model2>getMinimizationObjValue()>cutoff) { 

5696  // no good 

5697  delete model2; 

5698  delete solver; 

5699  return 2; 

5700  } 

5701  // get back solution 

5702  model.originalModel(model2,false); 

5703  delete model2; 

5704  int status; 

5705  if (model.getMinimizationObjValue()<cutoff&&model.bestSolution()) { 

5706  double objValue = model.getObjValue(); 

5707  const double * solution = model.bestSolution(); 

5708  setBestSolution(CBC_TREE_SOL,objValue,solution); 

5709  status = 0; 

5710  } else { 

5711  status=2; 

5712  } 

5713  if (model.status()) 

5714  status ++ ; // not finished search 

5715  delete solver; 

5716  return status; 

5717  } 

5718  // Set a pointer to a row cut which will be added instead of normal branching. 

5719  void 

5720  CbcModel::setNextRowCut(const OsiRowCut & cut) 

5721  { 

5722  nextRowCut_=new OsiRowCut(cut); 

5723  nextRowCut_>setEffectiveness(COIN_DBL_MAX); // mark so will always stay 

5724  } 

5725  /* Process root node and return a strengthened model 

5726  

5727  The method assumes that initialSolve() has been called to solve the 

5728  LP relaxation. It processes the root node and then returns a pointer 

5729  to the strengthened model (or NULL if infeasible) 

5730  */ 

5731  OsiSolverInterface * 

5732  CbcModel::strengthenedModel() 

5733  { 

5734  /* 

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

5736  */ 

5737  /* 

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

5739  CbcSimpleInteger object for each integer variable. 

5740  */ 

5741  findIntegers(false) ; 

5742  /* 

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

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

5745  */ 

5746  synchronizeModel() ; 

5747  

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

5749  continuousObjective_ = COIN_DBL_MAX ; 

5750  /* 

5751  Solve the relaxation. 

5752  

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

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

5755  continuous relaxation. 

5756  */ 

5757  bool feasible = resolve() ; 

5758  /* 

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

5760  continue with processing the root node. 

5761  */ 

5762  if (!feasible) 

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

5764  return NULL; } 

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

5766  originalContinuousObjective_ = solver_>getObjValue(); 

5767  

5768  /* 

5769  Begin setup to process a feasible root node. 

5770  */ 

5771  bestObjective_ = CoinMin(bestObjective_,1.0e50) ; 

5772  numberSolutions_ = 0 ; 

5773  numberHeuristicSolutions_ = 0 ; 

5774  // Everything is minimization 

5775  double cutoff=getCutoff() ; 

5776  double direction = solver_>getObjSense() ; 

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

5778  messageHandler()>message(CBC_CUTOFF_WARNING1, 

5779  messages()) 

5780  << cutoff << cutoff << CoinMessageEol ; 

5781  if (cutoff > bestObjective_) 

5782  cutoff = bestObjective_ ; 

5783  setCutoff(cutoff) ; 

5784  /* 

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

5786  */ 

5787  int numberColumns = getNumCols() ; 

5788  if (!currentSolution_) 

5789  currentSolution_ = new double[numberColumns] ; 

5790  testSolution_=currentSolution_; 

5791  /* 

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

5793  constraint system (aka the continuous system). 

5794  */ 

5795  continuousSolver_ = solver_>clone() ; 

5796  numberRowsAtContinuous_ = getNumRows() ; 

5797  /* 

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

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

5800  installed. 

5801  */ 

5802  analyzeObjective() ; 

5803  /* 

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

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

5806  that produced a given cut. 

5807  */ 

5808  int maximumWhich = 1000 ; 

5809  int * whichGenerator = new int[maximumWhich] ; 

5810  maximumNumberCuts_ = 0 ; 

5811  currentNumberCuts_ = 0 ; 

5812  delete [] addedCuts_ ; 

5813  addedCuts_ = NULL ; 

5814  /* 

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

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

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

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

5819  

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

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

5822  removed any slack cuts (numberOldActiveCuts and numberNewCuts have been 

5823  adjusted accordingly). 

5824  

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

5826  

5827  */ 

5828  int iCutGenerator; 

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

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

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

5832  } 

5833  OsiCuts cuts ; 

5834  int numberOldActiveCuts = 0 ; 

5835  int numberNewCuts = 0 ; 

5836  { int iObject ; 

5837  int preferredWay ; 

5838  int numberUnsatisfied = 0 ; 

5839  memcpy(currentSolution_,solver_>getColSolution(), 

5840  numberColumns*sizeof(double)) ; 

5841  

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

5843  { double infeasibility = 

5844  object_[iObject]>infeasibility(preferredWay) ; 

5845  if (infeasibility) numberUnsatisfied++ ; } 

5846  if (numberUnsatisfied) 

5847  { feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_, 

5848  NULL,numberOldActiveCuts,numberNewCuts, 

5849  maximumWhich, whichGenerator) ; } } 

5850  /* 

5851  We've taken the continuous relaxation as far as we can. 

5852  */ 

5853  

5854  OsiSolverInterface * newSolver=NULL; 

5855  if (feasible) { 

5856  // make copy of current solver 

5857  newSolver = solver_>clone(); 

5858  } 

5859  /* 

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

5861  */ 

5862  delete [] whichGenerator ; 

5863  delete [] walkback_ ; 

5864  walkback_ = NULL ; 

5865  delete [] addedCuts_ ; 

5866  addedCuts_ = NULL ; 

5867  if (continuousSolver_) 

5868  { delete continuousSolver_ ; 

5869  continuousSolver_ = NULL ; } 

5870  /* 

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

5872  */ 

5873  globalCuts_= OsiCuts() ; 

5874  

5875  return newSolver; 

5876  } 

5877  // Just update objectiveValue 

5878  void CbcModel::setBestObjectiveValue( double objectiveValue) 

5879  { 

5880  bestObjective_=objectiveValue; 

5881  } 

5882  double 

5883  CbcModel::getBestPossibleObjValue() const 

5884  { 

5885  return CoinMin(bestPossibleObjective_,bestObjective_) * solver_>getObjSense() ; 

5886  } 

5887  // Make given rows (L or G) into global cuts and remove from lp 

5888  void 

5889  CbcModel::makeGlobalCuts(int number,const int * which) 

5890  { 

5891  const double * rowLower = solver_>getRowLower(); 

5892  const double * rowUpper = solver_>getRowUpper(); 

5893  

5894  int numberRows = solver_>getNumRows(); 

5895  

5896  // Row copy 

5897  const double * elementByRow = solver_>getMatrixByRow()>getElements(); 

5898  const int * column = solver_>getMatrixByRow()>getIndices(); 

5899  const CoinBigIndex * rowStart = solver_>getMatrixByRow()>getVectorStarts(); 

5900  const int * rowLength = solver_>getMatrixByRow()>getVectorLengths(); 

5901  

5902  // Not all rows may be good so we need new array 

5903  int * whichDelete = new int[numberRows]; 

5904  int nDelete=0; 

5905  for (int i=0;i<number;i++) { 

5906  int iRow = which[i]; 

5907  if (iRow>=0&&iRow<numberRows) { 

5908  if (rowLower[iRow]<1.0e20rowUpper[iRow]>1.0e20) { 

5909  whichDelete[nDelete++]=iRow; 

5910  OsiRowCut thisCut; 

5911  thisCut.setLb(rowLower[iRow]); 

5912  thisCut.setUb(rowUpper[iRow]); 

5913  int start = rowStart[iRow]; 

5914  thisCut.setRow(rowLength[iRow],column+start,elementByRow+start); 

5915  globalCuts_.insert(thisCut) ; 

5916  } 

5917  } 

5918  } 

5919  if (nDelete) 

5920  solver_>deleteRows(nDelete,whichDelete); 

5921  delete [] whichDelete; 

5922  } 

5923  void 

5924  CbcModel::setNodeComparison(CbcCompareBase * compare) 

5925  { 

5926  delete nodeCompare_; 

5927  nodeCompare_ = compare>clone(); 

5928  } 

5929  void 

5930  CbcModel::setNodeComparison(CbcCompareBase & compare) 

5931  { 

5932  delete nodeCompare_; 

5933  nodeCompare_ = compare.clone(); 

5934  } 

5935  // Set the strategy. Clones 

5936  void 

5937  CbcModel::setStrategy(CbcStrategy & strategy) 

5938  { 

5939  delete strategy_; 

5940  strategy_ = strategy.clone(); 

5941  } 

5942  // Increases usedInSolution for nonzeros 

5943  void 

5944  CbcModel::incrementUsed(const double * solution) 

5945  { 

5946  // might as well mark all including continuous 

5947  int numberColumns = solver_>getNumCols(); 

5948  for (int i=0;i<numberColumns;i++) { 

5949  if (solution[i]) 

5950  usedInSolution_[i]++; 

5951  } 

5952  } 

5953  // Are there numerical difficulties (for initialSolve) ? 

5954  bool 

5955  CbcModel::isInitialSolveAbandoned() const 

5956  { 

5957  if (status_!=1) { 

5958  return false; 

5959  } else { 

5960  return solver_>isAbandoned(); 

5961  } 

5962  } 

5963  // Is optimality proven (for initialSolve) ? 

5964  bool 

5965  CbcModel::isInitialSolveProvenOptimal() const 

5966  { 

5967  if (status_!=1) { 

5968  return originalContinuousObjective_<1.0e50; 

