[1271]  1  /* $Id: CbcModel.cpp 1627 20110326 15:13:57Z tkr $ */ 

[2]  2  // Copyright (C) 2002, International Business Machines 

 3  // Corporation and others. All Rights Reserved. 

 4  #if defined(_MSC_VER) 

 5  // Turn off compiler warning about long names 

 6  # pragma warning(disable:4786) 

 7  #endif 

[311]  8  

[325]  9  #include "CbcConfig.h" 

[311]  10  

[2]  11  #include <string> 

 12  //#define CBC_DEBUG 1 

 13  //#define CHECK_CUT_COUNTS 

[1412]  14  //#define CHECK_NODE 

[2]  15  //#define CHECK_NODE_FULL 

[264]  16  //#define NODE_LOG 

[640]  17  //#define GLOBAL_CUTS_JUST_POINTERS 

[1271]  18  #ifdef CGL_DEBUG_GOMORY 

 19  extern int gomory_try; 

[687]  20  #endif 

[2]  21  #include <cassert> 

 22  #include <cmath> 

 23  #include <cfloat> 

[277]  24  

[311]  25  #ifdef COIN_HAS_CLP 

[208]  26  // include Presolve from Clp 

 27  #include "ClpPresolve.hpp" 

 28  #include "OsiClpSolverInterface.hpp" 

[931]  29  #include "ClpNode.hpp" 

 30  #include "ClpDualRowDantzig.hpp" 

 31  #include "ClpSimplexPrimal.hpp" 

[277]  32  #endif 

 33  

 34  #include "CbcEventHandler.hpp" 

[2]  35  

 36  #include "OsiSolverInterface.hpp" 

[264]  37  #include "OsiAuxInfo.hpp" 

[222]  38  #include "OsiSolverBranch.hpp" 

[640]  39  #include "OsiChooseVariable.hpp" 

[2]  40  #include "CoinWarmStartBasis.hpp" 

 41  #include "CoinPackedMatrix.hpp" 

[140]  42  #include "CoinHelperFunctions.hpp" 

[2]  43  #include "CbcBranchActual.hpp" 

[135]  44  #include "CbcBranchDynamic.hpp" 

[2]  45  #include "CbcHeuristic.hpp" 

[687]  46  #include "CbcHeuristicFPump.hpp" 

[1315]  47  #include "CbcHeuristicRINS.hpp" 

[1271]  48  #include "CbcHeuristicDive.hpp" 

[2]  49  #include "CbcModel.hpp" 

[186]  50  #include "CbcTreeLocal.hpp" 

[150]  51  #include "CbcStatistics.hpp" 

[97]  52  #include "CbcStrategy.hpp" 

[2]  53  #include "CbcMessage.hpp" 

 54  #include "OsiRowCut.hpp" 

 55  #include "OsiColCut.hpp" 

 56  #include "OsiRowCutDebugger.hpp" 

 57  #include "OsiCuts.hpp" 

 58  #include "CbcCountRowCut.hpp" 

 59  #include "CbcCutGenerator.hpp" 

[165]  60  #include "CbcFeasibilityBase.hpp" 

[687]  61  #include "CbcFathom.hpp" 

[2]  62  // include Probing 

 63  #include "CglProbing.hpp" 

[1271]  64  #include "CglGomory.hpp" 

 65  #include "CglTwomir.hpp" 

[222]  66  // include preprocessing 

 67  #include "CglPreProcess.hpp" 

[640]  68  #include "CglDuplicateRow.hpp" 

 69  #include "CglStored.hpp" 

[833]  70  #include "CglClique.hpp" 

[2]  71  

 72  #include "CoinTime.hpp" 

[640]  73  #include "CoinMpsIO.hpp" 

[2]  74  

 75  #include "CbcCompareActual.hpp" 

 76  #include "CbcTree.hpp" 

[1409]  77  // This may be dummy 

 78  #include "CbcThread.hpp" 

[2]  79  /* Various functions local to CbcModel.cpp */ 

 80  

 81  namespace { 

 82  

 83  // 

[1286]  84  // Returns the greatest common denominator of two 

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

[2]  86  // 

[1286]  87  static int gcd(int a, int b) 

[2]  88  { 

[1286]  89  int remainder = 1; 

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

 91  if (a > b) { 

 92  // Swap a and b 

 93  int temp = a; 

 94  a = b; 

 95  b = temp; 

[2]  96  } 

[1286]  97  // if zero then gcd is nonzero (zero may occur in rhs of packed) 

 98  if (!a) { 

 99  if (b) { 

 100  return b; 

 101  } else { 

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

 103  abort(); 

 104  } 

 105  } 

 106  while (remainder) { 

 107  remainder = b % a; 

 108  b = a; 

 109  a = remainder; 

 110  } 

 111  return b; 

[2]  112  } 

 113  

[931]  114  

[2]  115  

 116  #ifdef CHECK_NODE_FULL 

 117  

 118  /* 

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

 120  is 

 121  

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

 123  

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

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

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

 127  made to see if the numberPointingToThis_ field agrees. 

 128  */ 

 129  

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

 131  

[1286]  132  { 

 133  if (model.getNodeCount() == 661) return; 

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

 135  

 136  int j ; 

 137  int nNodes = branchingTree>size() ; 

[2]  138  # define MAXINFO 1000 

[1286]  139  int *count = new int [MAXINFO] ; 

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

 141  int nInfo = 0 ; 

 142  /* 

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

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

 145  branches remaining. 

[2]  146  

[1286]  147  TODO: The `while (nodeInfo)' loop could be made to break on reaching a 

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

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

 150  common ancestor. 

 151  */ 

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

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

 154  if (!node) 

 155  continue; 

 156  CbcNodeInfo *nodeInfo = node>nodeInfo() ; 

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

 158  assert(change) ; 

 159  while (nodeInfo) { 

 160  int k ; 

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

 162  if (nodeInfo == info[k]) break ; 

 163  } 

 164  if (k == nInfo) { 

 165  assert(nInfo < MAXINFO) ; 

 166  nInfo++ ; 

 167  info[k] = nodeInfo ; 

 168  count[k] = 0 ; 

 169  } 

 170  nodeInfo = nodeInfo>parent() ; 

 171  } 

 172  } 

 173  /* 

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

 175  increment the corresponding count. 

 176  */ 

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

 178  CbcNodeInfo *nodeInfo = info[j] ; 

 179  nodeInfo = nodeInfo>parent() ; 

 180  if (nodeInfo) { 

 181  int k ; 

 182  for (k = 0 ; k < nInfo ; k++) { 

 183  if (nodeInfo == info[k]) break ; 

 184  } 

 185  assert (k < nInfo) ; 

 186  count[k]++ ; 

 187  } 

 188  } 

 189  /* 

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

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

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

 193  references (unexplored branches, numberBranchesLeft()). 

 194  */ 

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

 196  CbcNodeInfo * nodeInfo = info[j] ; 

 197  if (nodeInfo) { 

 198  int k ; 

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

 200  if (nodeInfo == info[k]) 

 201  break ; 

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

 203  nodeInfo, 

 204  nodeInfo>numberBranchesLeft(), 

 205  nodeInfo>numberPointingToThis()) ; 

 206  assert(nodeInfo>numberPointingToThis() == 

 207  count[k] + nodeInfo>numberBranchesLeft()) ; 

 208  } 

 209  } 

[2]  210  

[1286]  211  delete [] count ; 

 212  delete [] info ; 

[2]  213  

[1286]  214  return ; 

 215  } 

 216  

[36]  217  #endif /* CHECK_NODE_FULL */ 

[2]  218  

[931]  219  

[2]  220  

 221  #ifdef CHECK_CUT_COUNTS 

 222  

 223  /* 

 224  Routine to verify that cut reference counts are correct. 

 225  */ 

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

 227  

[1286]  228  { 

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

[2]  230  

[1286]  231  int j ; 

 232  int nNodes = branchingTree>size() ; 

[2]  233  

[1286]  234  /* 

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

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

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

 238  */ 

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

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

 241  CbcNodeInfo * nodeInfo = node>nodeInfo() ; 

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

 243  while (nodeInfo) { 

 244  int k ; 

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

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

 247  if (cut) cut>tempNumber_ = 0; 

 248  } 

 249  nodeInfo = nodeInfo>parent() ; 

 250  } 

 251  } 

 252  /* 

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

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

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

 256  */ 

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

 258  CoinWarmStartBasis *debugws = model.getEmptyBasis() ; 

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

 260  CbcNodeInfo *nodeInfo = node>nodeInfo(); 

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

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

 263  node>nodeInfo(), node>columnNumber(), node>way(), 

 264  node>objectiveValue()) ; 

[2]  265  

[1286]  266  model.addCuts1(node, debugws) ; 

[2]  267  

[1286]  268  int i ; 

 269  int numberRowsAtContinuous = model.numberRowsAtContinuous() ; 

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

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

 272  CoinWarmStartBasis::Status status = 

 273  debugws>getArtifStatus(i + numberRowsAtContinuous) ; 

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

 275  addedCuts[i]>tempNumber_ += change ; 

 276  } 

 277  } 

[2]  278  

[1286]  279  while (nodeInfo) { 

 280  nodeInfo = nodeInfo>parent() ; 

 281  if (nodeInfo) printf(" > %x", nodeInfo); 

 282  } 

 283  printf("\n") ; 

 284  delete debugws ; 

 285  } 

 286  /* 

 287  The moment of truth: We've tallied up the references by direct scan of the search tree. Check for agreement with the count in the cut. 

[2]  288  

[1286]  289  TODO: Rewrite to check and print mismatch only when tempNumber_ == 0? 

 290  */ 

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

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

 293  CbcNodeInfo *nodeInfo = node>nodeInfo(); 

 294  while (nodeInfo) { 

 295  int k ; 

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

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

 298  if (cut && cut>tempNumber_ >= 0) { 

 299  if (cut>tempNumber_ != cut>numberPointingToThis()) 

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

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

 302  else 

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

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

 305  cut>tempNumber_ = 1 ; 

 306  } 

 307  } 

 308  nodeInfo = nodeInfo>parent() ; 

 309  } 

 310  } 

[2]  311  

[1286]  312  return ; 

 313  } 

[2]  314  

 315  #endif /* CHECK_CUT_COUNTS */ 

 316  

[931]  317  

[640]  318  #ifdef CHECK_CUT_SIZE 

 319  

 320  /* 

 321  Routine to verify that cut reference counts are correct. 

 322  */ 

 323  void verifyCutSize (const CbcTree * branchingTree, CbcModel &model) 

[1286]  324  { 

[640]  325  

[1286]  326  int j ; 

 327  int nNodes = branchingTree>size() ; 

 328  int totalCuts = 0; 

[640]  329  

[1286]  330  /* 

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

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

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

 334  */ 

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

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

 337  CbcNodeInfo * nodeInfo = node>nodeInfo() ; 

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

 339  while (nodeInfo) { 

 340  totalCuts += nodeInfo>numberCuts(); 

 341  nodeInfo = nodeInfo>parent() ; 

 342  } 

[640]  343  } 

[1286]  344  printf("*** CHECKING cuts (size) after %d nodes  %d cuts\n", model.getNodeCount(), totalCuts) ; 

 345  return ; 

[2]  346  } 

 347  

[640]  348  #endif /* CHECK_CUT_SIZE */ 

 349  

 350  } 

 351  

[1286]  352  /* End unnamed namespace for CbcModel.cpp */ 

[2]  353  

[931]  354  

[1286]  355  void 

[2]  356  CbcModel::analyzeObjective () 

 357  /* 

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

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

 360  and grabs the largest objective coefficient associated with an unfixed 

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

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

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

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

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

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

 367  

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

 369  

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

 371  */ 

[1286]  372  { 

 373  const double *objective = getObjCoefficients() ; 

 374  const double *lower = getColLower() ; 

 375  const double *upper = getColUpper() ; 

 376  /* 

 377  Scan continuous and integer variables to see if continuous 

 378  are cover or network with integral rhs. 

 379  */ 

 380  double continuousMultiplier = 1.0; 

 381  double * coeffMultiplier = NULL; 

 382  double largestObj = 0.0; 

 383  double smallestObj = COIN_DBL_MAX; 

 384  { 

 385  const double *rowLower = getRowLower() ; 

 386  const double *rowUpper = getRowUpper() ; 

 387  int numberRows = solver_>getNumRows() ; 

 388  double * rhs = new double [numberRows]; 

 389  memset(rhs, 0, numberRows*sizeof(double)); 

 390  int iColumn; 

 391  int numberColumns = solver_>getNumCols() ; 

 392  // Column copy of matrix 

 393  bool allPlusOnes = true; 

 394  bool allOnes = true; 

 395  int problemType = 1; 

 396  const double * element = solver_>getMatrixByCol()>getElements(); 

 397  const int * row = solver_>getMatrixByCol()>getIndices(); 

 398  const CoinBigIndex * columnStart = solver_>getMatrixByCol()>getVectorStarts(); 

 399  const int * columnLength = solver_>getMatrixByCol()>getVectorLengths(); 

 400  int numberInteger = 0; 

 401  int numberIntegerObj = 0; 

 402  int numberGeneralIntegerObj = 0; 

 403  int numberIntegerWeight = 0; 

 404  int numberContinuousObj = 0; 

 405  double cost = COIN_DBL_MAX; 

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

 407  if (upper[iColumn] == lower[iColumn]) { 

 408  CoinBigIndex start = columnStart[iColumn]; 

 409  CoinBigIndex end = start + columnLength[iColumn]; 

 410  for (CoinBigIndex j = start; j < end; j++) { 

 411  int iRow = row[j]; 

 412  rhs[iRow] += lower[iColumn] * element[j]; 

 413  } 

 414  } else { 

 415  double objValue = objective[iColumn]; 

 416  if (solver_>isInteger(iColumn)) 

 417  numberInteger++; 

 418  if (objValue) { 

 419  if (!solver_>isInteger(iColumn)) { 

 420  numberContinuousObj++; 

 421  } else { 

 422  largestObj = CoinMax(largestObj, fabs(objValue)); 

 423  smallestObj = CoinMin(smallestObj, fabs(objValue)); 

 424  numberIntegerObj++; 

 425  if (cost == COIN_DBL_MAX) 

 426  cost = objValue; 

 427  else if (cost != objValue) 

 428  cost = COIN_DBL_MAX; 

 429  int gap = static_cast<int> (upper[iColumn]  lower[iColumn]); 

 430  if (gap > 1) { 

 431  numberGeneralIntegerObj++; 

 432  numberIntegerWeight += gap; 

 433  } 

 434  } 

 435  } 

 436  } 

 437  } 

 438  int iType = 0; 

 439  if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 && 

 440  numberIntegerObj*3 < numberObjects_ && !parentModel_ && solver_>getNumRows() > 100) 

 441  iType = 3 + 4; 

 442  else if (!numberContinuousObj && numberIntegerObj <= 100 && 

 443  numberIntegerObj*5 < numberObjects_ && numberIntegerWeight <= 100 && 

 444  !parentModel_ && 

 445  solver_>getNumRows() > 100 && cost != COIN_DBL_MAX) 

 446  iType = 2 + 4; 

 447  else if (!numberContinuousObj && numberIntegerObj <= 100 && 

 448  numberIntegerObj*5 < numberObjects_ && 

 449  !parentModel_ && 

 450  solver_>getNumRows() > 100 && cost != COIN_DBL_MAX) 

 451  iType = 8; 

 452  int iTest = getMaximumNodes(); 

 453  if (iTest >= 987654320 && iTest < 987654330 && numberObjects_ && !parentModel_) { 

 454  iType = iTest  987654320; 

 455  printf("Testing %d integer variables out of %d objects (%d integer) have cost of %g  %d continuous\n", 

 456  numberIntegerObj, numberObjects_, numberInteger, cost, numberContinuousObj); 

 457  if (iType == 9) 

 458  exit(77); 

 459  if (numberContinuousObj) 

 460  iType = 0; 

 461  } 

 462  

 463  //if (!numberContinuousObj&&(numberIntegerObj<=5cost!=COIN_DBL_MAX)&& 

 464  //numberIntegerObj*3<numberObjects_&&!parentModel_&&solver_>getNumRows()>100) { 

 465  if (iType) { 

 466  /* 

 467  A) put high priority on (if none) 

 468  B) create artificial objective (if clp) 

 469  */ 

 470  int iPriority = 1; 

 471  for (int i = 0; i < numberObjects_; i++) { 

 472  int k = object_[i]>priority(); 

 473  if (iPriority == 1) 

 474  iPriority = k; 

 475  else if (iPriority != k) 

 476  iPriority = 2; 

 477  } 

 478  bool branchOnSatisfied = ((iType & 1) != 0); 

 479  bool createFake = ((iType & 2) != 0); 

 480  bool randomCost = ((iType & 4) != 0); 

 481  if (iPriority >= 0) { 

 482  char general[200]; 

 483  if (cost == COIN_DBL_MAX) { 

 484  sprintf(general, "%d integer variables out of %d objects (%d integer) have costs  high priority", 

 485  numberIntegerObj, numberObjects_, numberInteger); 

 486  } else if (cost == COIN_DBL_MAX) { 

 487  sprintf(general, "No integer variables out of %d objects (%d integer) have costs", 

 488  numberObjects_, numberInteger); 

 489  branchOnSatisfied = false; 

 490  } else { 

 491  sprintf(general, "%d integer variables out of %d objects (%d integer) have cost of %g  high priority", 

 492  numberIntegerObj, numberObjects_, numberInteger, cost); 

 493  } 

 494  messageHandler()>message(CBC_GENERAL, 

 495  messages()) 

 496  << general << CoinMessageEol ; 

 497  sprintf(general, "branch on satisfied %c create fake objective %c random cost %c", 

 498  branchOnSatisfied ? 'Y' : 'N', 

 499  createFake ? 'Y' : 'N', 

 500  randomCost ? 'Y' : 'N'); 

 501  messageHandler()>message(CBC_GENERAL, 

 502  messages()) 

 503  << general << CoinMessageEol ; 

 504  // switch off clp type branching 

 505  fastNodeDepth_ = 1; 

 506  int highPriority = (branchOnSatisfied) ? 999 : 100; 

 507  for (int i = 0; i < numberObjects_; i++) { 

 508  CbcSimpleInteger * thisOne = dynamic_cast <CbcSimpleInteger *> (object_[i]); 

 509  object_[i]>setPriority(1000); 

 510  if (thisOne) { 

 511  int iColumn = thisOne>columnNumber(); 

 512  if (objective[iColumn]) 

 513  thisOne>setPriority(highPriority); 

 514  } 

 515  } 

 516  } 

[1088]  517  #ifdef COIN_HAS_CLP 

[1286]  518  OsiClpSolverInterface * clpSolver 

 519  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 520  if (clpSolver && createFake) { 

 521  // Create artificial objective to be used when all else fixed 

 522  int numberColumns = clpSolver>getNumCols(); 

 523  double * fakeObj = new double [numberColumns]; 

 524  // Column copy 

 525  const CoinPackedMatrix * matrixByCol = clpSolver>getMatrixByCol(); 

 526  //const double * element = matrixByCol.getElements(); 

 527  //const int * row = matrixByCol.getIndices(); 

 528  //const CoinBigIndex * columnStart = matrixByCol.getVectorStarts(); 

 529  const int * columnLength = matrixByCol>getVectorLengths(); 

 530  const double * solution = clpSolver>getColSolution(); 

[1393]  531  #ifdef JJF_ZERO 

[1286]  532  int nAtBound = 0; 

 533  for (int i = 0; i < numberColumns; i++) { 

 534  double lowerValue = lower[i]; 

 535  double upperValue = upper[i]; 

 536  if (clpSolver>isInteger(i)) { 

 537  double lowerValue = lower[i]; 

 538  double upperValue = upper[i]; 

 539  double value = solution[i]; 

 540  if (value < lowerValue + 1.0e6  

 541  value > upperValue  1.0e6) 

 542  nAtBound++; 

 543  } 

 544  } 

[1271]  545  #endif 

[1409]  546  /* 

 547  Generate a random objective function for problems where the given objective 

 548  function is not terribly useful. (Nearly feasible, single integer variable, 

 549  that sort of thing. 

 550  */ 

[1286]  551  CoinDrand48(true, 1234567); 

 552  for (int i = 0; i < numberColumns; i++) { 

 553  double lowerValue = lower[i]; 

 554  double upperValue = upper[i]; 

 555  double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000) 

 556  : i + 1 + columnLength[i] * 1000; 

 557  value *= 0.001; 

 558  //value += columnLength[i]; 

 559  if (lowerValue > 1.0e5  upperValue < 1.0e5) { 

 560  if (fabs(lowerValue) > fabs(upperValue)) 

 561  value =  value; 

 562  if (clpSolver>isInteger(i)) { 

 563  double solValue = solution[i]; 

 564  // Better to add in 0.5 or 1.0?? 

 565  if (solValue < lowerValue + 1.0e6) 

 566  value = fabs(value) + 0.5; //fabs(value*1.5); 

 567  else if (solValue > upperValue  1.0e6) 

 568  value = fabs(value)  0.5; //fabs(value*1.5); 

 569  } 

 570  } else { 

 571  value = 0.0; 

 572  } 

 573  fakeObj[i] = value; 

 574  } 

 575  // pass to solver 

 576  clpSolver>setFakeObjective(fakeObj); 

 577  delete [] fakeObj; 

 578  } 

[1088]  579  #endif 

[1286]  580  } else if (largestObj < smallestObj*5.0 && !parentModel_ && 

 581  !numberContinuousObj && 

 582  !numberGeneralIntegerObj && 

 583  numberIntegerObj*2 < numberColumns) { 

 584  // up priorities on costed 

 585  int iPriority = 1; 

 586  for (int i = 0; i < numberObjects_; i++) { 

 587  int k = object_[i]>priority(); 

 588  if (iPriority == 1) 

 589  iPriority = k; 

 590  else if (iPriority != k) 

 591  iPriority = 2; 

 592  } 

 593  if (iPriority >= 100) { 

[1271]  594  #ifdef CLP_INVESTIGATE 

[1286]  595  printf("Setting variables with obj to high priority\n"); 

[1271]  596  #endif 

[1286]  597  for (int i = 0; i < numberObjects_; i++) { 

 598  CbcSimpleInteger * obj = 

 599  dynamic_cast <CbcSimpleInteger *>(object_[i]) ; 

 600  if (obj) { 

 601  int iColumn = obj>columnNumber(); 

 602  if (objective[iColumn]) 

 603  object_[i]>setPriority(iPriority  1); 

 604  } 

 605  } 

 606  } 

 607  } 

 608  int iRow; 

 609  for (iRow = 0; iRow < numberRows; iRow++) { 

 610  if (rowLower[iRow] > 1.0e20 && 

 611  fabs(rowLower[iRow]  rhs[iRow]  floor(rowLower[iRow]  rhs[iRow] + 0.5)) > 1.0e10) { 

 612  continuousMultiplier = 0.0; 

 613  break; 

 614  } 

 615  if (rowUpper[iRow] < 1.0e20 && 

 616  fabs(rowUpper[iRow]  rhs[iRow]  floor(rowUpper[iRow]  rhs[iRow] + 0.5)) > 1.0e10) { 

 617  continuousMultiplier = 0.0; 

 618  break; 

 619  } 

 620  // set rhs to limiting value 

 621  if (rowLower[iRow] != rowUpper[iRow]) { 

 622  if (rowLower[iRow] > 1.0e20) { 

 623  if (rowUpper[iRow] < 1.0e20) { 

 624  // no good 

 625  continuousMultiplier = 0.0; 

 626  break; 

 627  } else { 

 628  rhs[iRow] = rowLower[iRow]  rhs[iRow]; 

 629  if (problemType < 0) 

 630  problemType = 3; // set cover 

 631  else if (problemType != 3) 

 632  problemType = 4; 

 633  } 

 634  } else { 

 635  rhs[iRow] = rowUpper[iRow]  rhs[iRow]; 

 636  if (problemType < 0) 

 637  problemType = 1; // set partitioning <= 

 638  else if (problemType != 1) 

 639  problemType = 4; 

 640  } 

 641  } else { 

 642  rhs[iRow] = rowUpper[iRow]  rhs[iRow]; 

 643  if (problemType < 0) 

 644  problemType = 3; // set partitioning == 

 645  else if (problemType != 2) 

 646  problemType = 2; 

 647  } 

 648  if (fabs(rhs[iRow]  1.0) > 1.0e12) 

 649  problemType = 4; 

 650  } 

 651  if (continuousMultiplier) { 

 652  // 1 network, 2 cover, 4 negative cover 

 653  int possible = 7; 

 654  bool unitRhs = true; 

 655  // See which rows could be set cover 

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

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

 658  CoinBigIndex start = columnStart[iColumn]; 

 659  CoinBigIndex end = start + columnLength[iColumn]; 

 660  for (CoinBigIndex j = start; j < end; j++) { 

 661  double value = element[j]; 

 662  if (value == 1.0) { 

 663  } else if (value == 1.0) { 

 664  rhs[row[j]] = 0.5; 

 665  allPlusOnes = false; 

 666  } else { 

 667  rhs[row[j]] = COIN_DBL_MAX; 

 668  allOnes = false; 

 669  } 

 670  } 

 671  } 

 672  } 

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

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

 675  if (!isInteger(iColumn)) { 

 676  CoinBigIndex start = columnStart[iColumn]; 

 677  CoinBigIndex end = start + columnLength[iColumn]; 

 678  double rhsValue = 0.0; 

 679  // 1 all ones, 1 all 1s, 2 all + 1, 3 no good 

 680  int type = 0; 

 681  for (CoinBigIndex j = start; j < end; j++) { 

 682  double value = element[j]; 

 683  if (fabs(value) != 1.0) { 

 684  type = 3; 

 685  break; 

 686  } else if (value == 1.0) { 

 687  if (!type) 

 688  type = 1; 

 689  else if (type != 1) 

 690  type = 2; 

 691  } else { 

 692  if (!type) 

 693  type = 1; 

 694  else if (type != 1) 

 695  type = 2; 

 696  } 

 697  int iRow = row[j]; 

 698  if (rhs[iRow] == COIN_DBL_MAX) { 

 699  type = 3; 

 700  break; 

 701  } else if (rhs[iRow] == 0.5) { 

 702  // different values 

 703  unitRhs = false; 

 704  } else if (rhsValue) { 

 705  if (rhsValue != rhs[iRow]) 

 706  unitRhs = false; 

 707  } else { 

 708  rhsValue = rhs[iRow]; 

 709  } 

 710  } 

 711  // if no elements OK 

 712  if (type == 3) { 

 713  // no good 

 714  possible = 0; 

 715  break; 

 716  } else if (type == 2) { 

 717  if (end  start > 2) { 

 718  // no good 

 719  possible = 0; 

 720  break; 

 721  } else { 

 722  // only network 

 723  possible &= 1; 

 724  if (!possible) 

 725  break; 

 726  } 

 727  } else if (type == 1) { 

 728  // only cover 

 729  possible &= 2; 

 730  if (!possible) 

 731  break; 

 732  } else if (type == 1) { 

 733  // only negative cover 

 734  possible &= 4; 

 735  if (!possible) 

 736  break; 

 737  } 

 738  } 

 739  } 

 740  } 

 741  if ((possible == 2  possible == 4) && !unitRhs) { 

[931]  742  #if COIN_DEVELOP>1 

[1286]  743  printf("XXXXXX Continuous all +1 but different rhs\n"); 

[810]  744  #endif 

[1286]  745  possible = 0; 

 746  } 

 747  // may be all integer 

 748  if (possible != 7) { 

 749  if (!possible) 

 750  continuousMultiplier = 0.0; 

 751  else if (possible == 1) 

 752  continuousMultiplier = 1.0; 

 753  else 

 754  continuousMultiplier = 0.0; // 0.5 was incorrect; 

[931]  755  #if COIN_DEVELOP>1 

[1286]  756  if (continuousMultiplier) 

 757  printf("XXXXXX multiplier of %g\n", continuousMultiplier); 

[810]  758  #endif 

[1286]  759  if (continuousMultiplier == 0.5) { 

 760  coeffMultiplier = new double [numberColumns]; 

 761  bool allOne = true; 

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

 763  coeffMultiplier[iColumn] = 1.0; 

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

 765  if (!isInteger(iColumn)) { 

 766  CoinBigIndex start = columnStart[iColumn]; 

 767  int iRow = row[start]; 

 768  double value = rhs[iRow]; 

 769  assert (value >= 0.0); 

 770  if (value != 0.0 && value != 1.0) 

 771  allOne = false; 

 772  coeffMultiplier[iColumn] = 0.5 * value; 

 773  } 

 774  } 

 775  } 

 776  if (allOne) { 

 777  // back to old way 

 778  delete [] coeffMultiplier; 

 779  coeffMultiplier = NULL; 

 780  } 

 781  } 

 782  } else { 

 783  // all integer 

 784  problemType_ = problemType; 

[931]  785  #if COIN_DEVELOP>1 

[1286]  786  printf("Problem type is %d\n", problemType_); 

[814]  787  #endif 

[1286]  788  } 

 789  } 

 790  

 791  // But try again 

 792  if (continuousMultiplier < 1.0) { 

 793  memset(rhs, 0, numberRows*sizeof(double)); 

 794  int * count = new int [numberRows]; 

 795  memset(count, 0, numberRows*sizeof(int)); 

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

 797  CoinBigIndex start = columnStart[iColumn]; 

 798  CoinBigIndex end = start + columnLength[iColumn]; 

 799  if (upper[iColumn] == lower[iColumn]) { 

 800  for (CoinBigIndex j = start; j < end; j++) { 

 801  int iRow = row[j]; 

 802  rhs[iRow] += lower[iColumn] * element[j]; 

 803  } 

 804  } else if (solver_>isInteger(iColumn)) { 

 805  for (CoinBigIndex j = start; j < end; j++) { 

 806  int iRow = row[j]; 

 807  if (fabs(element[j]  floor(element[j] + 0.5)) > 1.0e10) 

 808  rhs[iRow] = COIN_DBL_MAX; 

 809  } 

 810  } else { 

 811  for (CoinBigIndex j = start; j < end; j++) { 

 812  int iRow = row[j]; 

 813  count[iRow]++; 

 814  if (fabs(element[j]) != 1.0) 

 815  rhs[iRow] = COIN_DBL_MAX; 

 816  } 

 817  } 

 818  } 

 819  // now look at continuous 

 820  bool allGood = true; 

 821  double direction = solver_>getObjSense() ; 

 822  int numberObj = 0; 

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

 824  if (upper[iColumn] > lower[iColumn]) { 

 825  double objValue = objective[iColumn] * direction; 

 826  if (objValue && !solver_>isInteger(iColumn)) { 

 827  numberObj++; 

 828  CoinBigIndex start = columnStart[iColumn]; 

 829  CoinBigIndex end = start + columnLength[iColumn]; 

 830  if (objValue > 0.0) { 

 831  // wants to be as low as possible 

 832  if (lower[iColumn] < 1.0e10  fabs(lower[iColumn]  floor(lower[iColumn] + 0.5)) > 1.0e10) { 

 833  allGood = false; 

 834  break; 

 835  } else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn]  floor(upper[iColumn] + 0.5)) > 1.0e10) { 

 836  allGood = false; 

 837  break; 

 838  } 

 839  bool singletonRow = true; 

 840  bool equality = false; 

 841  for (CoinBigIndex j = start; j < end; j++) { 

 842  int iRow = row[j]; 

 843  if (count[iRow] > 1) 

 844  singletonRow = false; 

 845  else if (rowLower[iRow] == rowUpper[iRow]) 

 846  equality = true; 

 847  double rhsValue = rhs[iRow]; 

 848  double lowerValue = rowLower[iRow]; 

 849  double upperValue = rowUpper[iRow]; 

 850  if (rhsValue < 1.0e20) { 

 851  if (lowerValue > 1.0e20) 

 852  lowerValue = rhsValue; 

 853  if (upperValue < 1.0e20) 

 854  upperValue = rhsValue; 

 855  } 

 856  if (fabs(rhsValue) > 1.0e20  fabs(rhsValue  floor(rhsValue + 0.5)) > 1.0e10 

 857   fabs(element[j]) != 1.0) { 

 858  // no good 

 859  allGood = false; 

 860  break; 

 861  } 

 862  if (element[j] > 0.0) { 

 863  if (lowerValue > 1.0e20 && fabs(lowerValue  floor(lowerValue + 0.5)) > 1.0e10) { 

 864  // no good 

 865  allGood = false; 

 866  break; 

 867  } 

 868  } else { 

 869  if (upperValue < 1.0e20 && fabs(upperValue  floor(upperValue + 0.5)) > 1.0e10) { 

 870  // no good 

 871  allGood = false; 

 872  break; 

 873  } 

 874  } 

 875  } 

 876  if (!singletonRow && end > start + 1 && !equality) 

 877  allGood = false; 

 878  if (!allGood) 

 879  break; 

 880  } else { 

 881  // wants to be as high as possible 

 882  if (upper[iColumn] > 1.0e10  fabs(upper[iColumn]  floor(upper[iColumn] + 0.5)) > 1.0e10) { 

 883  allGood = false; 

 884  break; 

 885  } else if (lower[iColumn] > 1.0e10 && fabs(lower[iColumn]  floor(lower[iColumn] + 0.5)) > 1.0e10) { 

 886  allGood = false; 

 887  break; 

 888  } 

 889  bool singletonRow = true; 

 890  bool equality = false; 

 891  for (CoinBigIndex j = start; j < end; j++) { 

 892  int iRow = row[j]; 

 893  if (count[iRow] > 1) 

 894  singletonRow = false; 

 895  else if (rowLower[iRow] == rowUpper[iRow]) 

 896  equality = true; 

 897  double rhsValue = rhs[iRow]; 

 898  double lowerValue = rowLower[iRow]; 

 899  double upperValue = rowUpper[iRow]; 

 900  if (rhsValue < 1.0e20) { 

 901  if (lowerValue > 1.0e20) 

 902  lowerValue = rhsValue; 

 903  if (upperValue < 1.0e20) 

 904  upperValue = rhsValue; 

 905  } 

 906  if (fabs(rhsValue) > 1.0e20  fabs(rhsValue  floor(rhsValue + 0.5)) > 1.0e10 

 907   fabs(element[j]) != 1.0) { 

 908  // no good 

 909  allGood = false; 

 910  break; 

 911  } 

 912  if (element[j] < 0.0) { 

 913  if (lowerValue > 1.0e20 && fabs(lowerValue  floor(lowerValue + 0.5)) > 1.0e10) { 

 914  // no good 

 915  allGood = false; 

 916  break; 

 917  } 

 918  } else { 

 919  if (upperValue < 1.0e20 && fabs(upperValue  floor(upperValue + 0.5)) > 1.0e10) { 

 920  // no good 

 921  allGood = false; 

 922  break; 

 923  } 

 924  } 

 925  } 

 926  if (!singletonRow && end > start + 1 && !equality) 

 927  allGood = false; 

 928  if (!allGood) 

 929  break; 

 930  } 

 931  } 

 932  } 

 933  } 

 934  delete [] count; 

 935  if (allGood) { 

[931]  936  #if COIN_DEVELOP>1 

[1286]  937  if (numberObj) 

 938  printf("YYYY analysis says all continuous with costs will be integer\n"); 

[837]  939  #endif 

[1286]  940  continuousMultiplier = 1.0; 

 941  } 

 942  } 

 943  delete [] rhs; 

[833]  944  } 

[1286]  945  /* 

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

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

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

[2]  949  

[1286]  950  If the objective is found to contain only integer variables, set the 

 951  fathoming discipline to strict. 

 952  */ 

 953  double maximumCost = 0.0 ; 

 954  //double trueIncrement=0.0; 

 955  int iColumn ; 

 956  int numberColumns = getNumCols() ; 

 957  double scaleFactor = 1.0; // due to rhs etc 

[1409]  958  /* 

 959  Original model did not have integer bounds. 

 960  */ 

[1286]  961  if ((specialOptions_&65536) == 0) { 

 962  /* be on safe side (later look carefully as may be able to 

 963  to get 0.5 say if bounds are multiples of 0.5 */ 

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

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

 966  double value; 

 967  value = fabs(lower[iColumn]); 

 968  if (floor(value + 0.5) != value) { 

 969  scaleFactor = CoinMin(scaleFactor, 0.5); 

 970  if (floor(2.0*value + 0.5) != 2.0*value) { 

 971  scaleFactor = CoinMin(scaleFactor, 0.25); 

 972  if (floor(4.0*value + 0.5) != 4.0*value) { 

 973  scaleFactor = 0.0; 

 974  } 

 975  } 

 976  } 

 977  value = fabs(upper[iColumn]); 

 978  if (floor(value + 0.5) != value) { 

 979  scaleFactor = CoinMin(scaleFactor, 0.5); 

 980  if (floor(2.0*value + 0.5) != 2.0*value) { 

 981  scaleFactor = CoinMin(scaleFactor, 0.25); 

 982  if (floor(4.0*value + 0.5) != 4.0*value) { 

 983  scaleFactor = 0.0; 

 984  } 

 985  } 

 986  } 

 987  } 

 988  } 

 989  } 

 990  bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0 ; 

 991  if (possibleMultiple) { 

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

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

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

 995  } 

 996  } 

 997  } 

 998  setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple) ; 

 999  /* 

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

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

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

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

 1004  later). 

[2]  1005  

[1286]  1006  2520.0 is used as it is a nice multiple of 2,3,5,7 

 1007  */ 

 1008  if (possibleMultiple && maximumCost) { 

 1009  int increment = 0 ; 

 1010  double multiplier = 2520.0 ; 

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

 1012  multiplier *= 10.0 ; 

 1013  int bigIntegers = 0; // Count of large costs which are integer 

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

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

 1016  double objValue = fabs(objective[iColumn]); 

 1017  if (!isInteger(iColumn)) { 

 1018  if (!coeffMultiplier) 

 1019  objValue *= continuousMultiplier; 

 1020  else 

 1021  objValue *= coeffMultiplier[iColumn]; 

 1022  } 

 1023  if (objValue) { 

 1024  double value = objValue * multiplier ; 

 1025  if (value < 2.1e9) { 

 1026  int nearest = static_cast<int> (floor(value + 0.5)) ; 

 1027  if (fabs(value  floor(value + 0.5)) > 1.0e8) { 

 1028  increment = 0 ; 

 1029  break ; 

 1030  } else if (!increment) { 

 1031  increment = nearest ; 

 1032  } else { 

 1033  increment = gcd(increment, nearest) ; 

 1034  } 

 1035  } else { 

 1036  // large value  may still be multiple of 1.0 

 1037  if (fabs(objValue  floor(objValue + 0.5)) > 1.0e8) { 

 1038  increment = 0; 

 1039  break; 

 1040  } else { 

 1041  bigIntegers++; 

 1042  } 

 1043  } 

 1044  } 

 1045  } 

 1046  } 

 1047  delete [] coeffMultiplier; 

 1048  /* 

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

 1050  */ 

 1051  if (increment) { 

 1052  double value = increment ; 

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

 1054  if (bigIntegers) { 

 1055  // allow for 1.0 

 1056  increment = gcd(increment, static_cast<int> (multiplier)); 

 1057  value = increment; 

 1058  } 

 1059  value /= multiplier ; 

 1060  value *= scaleFactor; 

 1061  //trueIncrement=CoinMax(cutoff,value);; 

 1062  if (value*0.999 > cutoff) { 

 1063  messageHandler()>message(CBC_INTEGERINCREMENT, 

 1064  messages()) 

 1065  << value << CoinMessageEol ; 

 1066  setDblParam(CbcModel::CbcCutoffIncrement, value*0.999) ; 

 1067  } 

 1068  } 

[474]  1069  } 

[2]  1070  

[1286]  1071  return ; 

[2]  1072  } 

 1073  

[1332]  1074  /* 

 1075  saveModel called (carved out of) BranchandBound 

 1076  */ 

[1330]  1077  void CbcModel::saveModel(OsiSolverInterface * saveSolver, double * checkCutoffForRestart, bool * feasible) 

 1078  { 

 1079  if (saveSolver && (specialOptions_&32768) != 0) { 

 1080  // See if worth trying reduction 

 1081  *checkCutoffForRestart = getCutoff(); 

 1082  bool tryNewSearch = solverCharacteristics_>reducedCostsAccurate() && 

 1083  (*checkCutoffForRestart < 1.0e20); 

 1084  int numberColumns = getNumCols(); 

 1085  if (tryNewSearch) { 

 1086  #ifdef CLP_INVESTIGATE 

 1087  printf("after %d nodes, cutoff %g  looking\n", 

 1088  numberNodes_, getCutoff()); 

 1089  #endif 

 1090  saveSolver>resolve(); 

 1091  double direction = saveSolver>getObjSense() ; 

 1092  double gap = *checkCutoffForRestart  saveSolver>getObjValue() * direction ; 

 1093  double tolerance; 

 1094  saveSolver>getDblParam(OsiDualTolerance, tolerance) ; 

 1095  if (gap <= 0.0) 

 1096  gap = tolerance; 

 1097  gap += 100.0 * tolerance; 

 1098  double integerTolerance = getDblParam(CbcIntegerTolerance) ; 

 1099  

 1100  const double *lower = saveSolver>getColLower() ; 

 1101  const double *upper = saveSolver>getColUpper() ; 

 1102  const double *solution = saveSolver>getColSolution() ; 

 1103  const double *reducedCost = saveSolver>getReducedCost() ; 

 1104  

 1105  int numberFixed = 0 ; 

 1106  int numberFixed2 = 0; 

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

 1108  int iColumn = integerVariable_[i] ; 

 1109  double djValue = direction * reducedCost[iColumn] ; 

 1110  if (upper[iColumn]  lower[iColumn] > integerTolerance) { 

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

 1112  saveSolver>setColUpper(iColumn, lower[iColumn]) ; 

 1113  numberFixed++ ; 

 1114  } else if (solution[iColumn] > upper[iColumn]  integerTolerance && djValue > gap) { 

 1115  saveSolver>setColLower(iColumn, upper[iColumn]) ; 

 1116  numberFixed++ ; 

 1117  } 

 1118  } else { 

 1119  numberFixed2++; 

 1120  } 

 1121  } 

 1122  #ifdef COIN_DEVELOP 

[1409]  1123  /* 

 1124  We're debugging. (specialOptions 1) 

 1125  */ 

[1330]  1126  if ((specialOptions_&1) != 0) { 

 1127  const OsiRowCutDebugger *debugger = saveSolver>getRowCutDebugger() ; 

 1128  if (debugger) { 

 1129  printf("Contains optimal\n") ; 

 1130  OsiSolverInterface * temp = saveSolver>clone(); 

 1131  const double * solution = debugger>optimalSolution(); 

 1132  const double *lower = temp>getColLower() ; 

 1133  const double *upper = temp>getColUpper() ; 

 1134  int n = temp>getNumCols(); 

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

 1136  if (temp>isInteger(i)) { 

 1137  double value = floor(solution[i] + 0.5); 

 1138  assert (value >= lower[i] && value <= upper[i]); 

 1139  temp>setColLower(i, value); 

 1140  temp>setColUpper(i, value); 

 1141  } 

 1142  } 

 1143  temp>writeMps("reduced_fix"); 

 1144  delete temp; 

 1145  saveSolver>writeMps("reduced"); 

 1146  } else { 

 1147  abort(); 

 1148  } 

 1149  } 

 1150  printf("Restart could fix %d integers (%d already fixed)\n", 

 1151  numberFixed + numberFixed2, numberFixed2); 

 1152  #endif 

 1153  numberFixed += numberFixed2; 

 1154  if (numberFixed*20 < numberColumns) 

 1155  tryNewSearch = false; 

 1156  } 

 1157  if (tryNewSearch) { 

 1158  // back to solver without cuts? 

 1159  OsiSolverInterface * solver2 = continuousSolver_>clone(); 

 1160  const double *lower = saveSolver>getColLower() ; 

 1161  const double *upper = saveSolver>getColUpper() ; 

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

 1163  int iColumn = integerVariable_[i] ; 

 1164  solver2>setColLower(iColumn, lower[iColumn]); 

 1165  solver2>setColUpper(iColumn, upper[iColumn]); 

 1166  } 

 1167  // swap 

 1168  delete saveSolver; 

 1169  saveSolver = solver2; 

 1170  double * newSolution = new double[numberColumns]; 

 1171  double objectiveValue = *checkCutoffForRestart; 

 1172  CbcSerendipity heuristic(*this); 

 1173  if (bestSolution_) 

 1174  heuristic.setInputSolution(bestSolution_, bestObjective_); 

 1175  heuristic.setFractionSmall(0.9); 

 1176  heuristic.setFeasibilityPumpOptions(1008013); 

 1177  // Use numberNodes to say how many are original rows 

 1178  heuristic.setNumberNodes(continuousSolver_>getNumRows()); 

 1179  #ifdef COIN_DEVELOP 

 1180  if (continuousSolver_>getNumRows() < 

 1181  saveSolver>getNumRows()) 

 1182  printf("%d rows added ZZZZZ\n", 

 1183  solver_>getNumRows()  continuousSolver_>getNumRows()); 

 1184  #endif 

 1185  int returnCode = heuristic.smallBranchAndBound(saveSolver, 

 1186  1, newSolution, 

 1187  objectiveValue, 

 1188  *checkCutoffForRestart, "Reduce"); 

 1189  if (returnCode < 0) { 

 1190  #ifdef COIN_DEVELOP 

 1191  printf("Restart  not small enough to do search after fixing\n"); 

 1192  #endif 

 1193  delete [] newSolution; 

 1194  } else { 

 1195  if ((returnCode&1) != 0) { 

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

 1197  numberSolutions_++; 

 1198  numberHeuristicSolutions_++; 

 1199  lastHeuristic_ = NULL; 

 1200  setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ; 

 1201  } 

 1202  delete [] newSolution; 

 1203  *feasible = false; // stop search 

 1204  } 

[1412]  1205  #if 0 // probably not needed def CBC_THREAD 

 1206  if (master_) { 

 1207  lockThread(); 

 1208  if (parallelMode() > 0) { 

 1209  while (master_>waitForThreadsInTree(0)) { 

 1210  lockThread(); 

 1211  double dummyBest; 

 1212  tree_>cleanTree(this, COIN_DBL_MAX, dummyBest) ; 

 1213  //unlockThread(); 

 1214  } 

 1215  } 

 1216  master_>waitForThreadsInTree(2); 

 1217  delete master_; 

 1218  master_ = NULL; 

 1219  masterThread_ = NULL; 

 1220  } 

 1221  #endif 

[1330]  1222  } 

 1223  } 

 1224  } 

 1225  /* 

 1226  Adds integers, called from BranchandBound() 

 1227  */ 

 1228  void CbcModel::AddIntegers() 

 1229  { 

 1230  int numberColumns = continuousSolver_>getNumCols(); 

 1231  int numberRows = continuousSolver_>getNumRows(); 

 1232  int * del = new int [CoinMax(numberColumns, numberRows)]; 

 1233  int * original = new int [numberColumns]; 

 1234  char * possibleRow = new char [numberRows]; 

 1235  { 

 1236  const CoinPackedMatrix * rowCopy = continuousSolver_>getMatrixByRow(); 

 1237  const int * column = rowCopy>getIndices(); 

 1238  const int * rowLength = rowCopy>getVectorLengths(); 

 1239  const CoinBigIndex * rowStart = rowCopy>getVectorStarts(); 

 1240  const double * rowLower = continuousSolver_>getRowLower(); 

 1241  const double * rowUpper = continuousSolver_>getRowUpper(); 

 1242  const double * element = rowCopy>getElements(); 

 1243  for (int i = 0; i < numberRows; i++) { 

 1244  int nLeft = 0; 

 1245  bool possible = false; 

 1246  if (rowLower[i] < 1.0e20) { 

 1247  double value = rowUpper[i]; 

 1248  if (fabs(value  floor(value + 0.5)) < 1.0e8) 

 1249  possible = true; 

 1250  } else if (rowUpper[i] > 1.0e20) { 

 1251  double value = rowLower[i]; 

 1252  if (fabs(value  floor(value + 0.5)) < 1.0e8) 

 1253  possible = true; 

 1254  } else { 

 1255  double value = rowUpper[i]; 

 1256  if (rowLower[i] == rowUpper[i] && 

 1257  fabs(value  floor(value + 0.5)) < 1.0e8) 

 1258  possible = true; 

 1259  } 

 1260  for (CoinBigIndex j = rowStart[i]; 

 1261  j < rowStart[i] + rowLength[i]; j++) { 

 1262  int iColumn = column[j]; 

 1263  if (continuousSolver_>isInteger(iColumn)) { 

 1264  if (fabs(element[j]) != 1.0) 

 1265  possible = false; 

 1266  } else { 

 1267  nLeft++; 

 1268  } 

 1269  } 

 1270  if (possible  !nLeft) 

 1271  possibleRow[i] = 1; 

 1272  else 

 1273  possibleRow[i] = 0; 

 1274  } 

 1275  } 

 1276  int nDel = 0; 

 1277  for (int i = 0; i < numberColumns; i++) { 

 1278  original[i] = i; 

 1279  if (continuousSolver_>isInteger(i)) 

 1280  del[nDel++] = i; 

 1281  } 

 1282  int nExtra = 0; 

 1283  OsiSolverInterface * copy1 = continuousSolver_>clone(); 

 1284  int nPass = 0; 

 1285  while (nDel && nPass < 10) { 

 1286  nPass++; 

 1287  OsiSolverInterface * copy2 = copy1>clone(); 

 1288  int nLeft = 0; 

 1289  for (int i = 0; i < nDel; i++) 

 1290  original[del[i]] = 1; 

 1291  for (int i = 0; i < numberColumns; i++) { 

 1292  int kOrig = original[i]; 

 1293  if (kOrig >= 0) 

 1294  original[nLeft++] = kOrig; 

 1295  } 

 1296  assert (nLeft == numberColumns  nDel); 

 1297  copy2>deleteCols(nDel, del); 

 1298  numberColumns = copy2>getNumCols(); 

 1299  const CoinPackedMatrix * rowCopy = copy2>getMatrixByRow(); 

 1300  numberRows = rowCopy>getNumRows(); 

 1301  const int * column = rowCopy>getIndices(); 

 1302  const int * rowLength = rowCopy>getVectorLengths(); 

 1303  const CoinBigIndex * rowStart = rowCopy>getVectorStarts(); 

 1304  const double * rowLower = copy2>getRowLower(); 

 1305  const double * rowUpper = copy2>getRowUpper(); 

 1306  const double * element = rowCopy>getElements(); 

 1307  const CoinPackedMatrix * columnCopy = copy2>getMatrixByCol(); 

 1308  const int * columnLength = columnCopy>getVectorLengths(); 

 1309  nDel = 0; 

 1310  // Could do gcd stuff on ones with costs 

 1311  for (int i = 0; i < numberRows; i++) { 

 1312  if (!rowLength[i]) { 

 1313  del[nDel++] = i; 

 1314  possibleRow[i] = 1; 

 1315  } else if (possibleRow[i]) { 

 1316  if (rowLength[i] == 1) { 

 1317  int k = rowStart[i]; 

 1318  int iColumn = column[k]; 

 1319  if (!copy2>isInteger(iColumn)) { 

 1320  double mult = 1.0 / fabs(element[k]); 

 1321  if (rowLower[i] < 1.0e20) { 

 1322  double value = rowUpper[i] * mult; 

 1323  if (fabs(value  floor(value + 0.5)) < 1.0e8) { 

 1324  del[nDel++] = i; 

 1325  if (columnLength[iColumn] == 1) { 

 1326  copy2>setInteger(iColumn); 

 1327  int kOrig = original[iColumn]; 

 1328  setOptionalInteger(kOrig); 

 1329  } 

 1330  } 

 1331  } else if (rowUpper[i] > 1.0e20) { 

 1332  double value = rowLower[i] * mult; 

 1333  if (fabs(value  floor(value + 0.5)) < 1.0e8) { 

 1334  del[nDel++] = i; 

 1335  if (columnLength[iColumn] == 1) { 

 1336  copy2>setInteger(iColumn); 

 1337  int kOrig = original[iColumn]; 

 1338  setOptionalInteger(kOrig); 

 1339  } 

 1340  } 

 1341  } else { 

 1342  double value = rowUpper[i] * mult; 

 1343  if (rowLower[i] == rowUpper[i] && 

 1344  fabs(value  floor(value + 0.5)) < 1.0e8) { 

 1345  del[nDel++] = i; 

 1346  copy2>setInteger(iColumn); 

 1347  int kOrig = original[iColumn]; 

 1348  setOptionalInteger(kOrig); 

 1349  } 

 1350  } 

 1351  } 

 1352  } else { 

 1353  // only if all singletons 

 1354  bool possible = false; 

 1355  if (rowLower[i] < 1.0e20) { 

 1356  double value = rowUpper[i]; 

 1357  if (fabs(value  floor(value + 0.5)) < 1.0e8) 

 1358  possible = true; 

 1359  } else if (rowUpper[i] > 1.0e20) { 

 1360  double value = rowLower[i]; 

 1361  if (fabs(value  floor(value + 0.5)) < 1.0e8) 

 1362  possible = true; 

 1363  } else { 

 1364  double value = rowUpper[i]; 

 1365  if (rowLower[i] == rowUpper[i] && 

 1366  fabs(value  floor(value + 0.5)) < 1.0e8) 

 1367  possible = true; 

 1368  } 

 1369  if (possible) { 

 1370  for (CoinBigIndex j = rowStart[i]; 

 1371  j < rowStart[i] + rowLength[i]; j++) { 

 1372  int iColumn = column[j]; 

 1373  if (columnLength[iColumn] != 1  fabs(element[j]) != 1.0) { 

 1374  possible = false; 

 1375  break; 

 1376  } 

 1377  } 

 1378  if (possible) { 

 1379  for (CoinBigIndex j = rowStart[i]; 

 1380  j < rowStart[i] + rowLength[i]; j++) { 

 1381  int iColumn = column[j]; 

 1382  if (!copy2>isInteger(iColumn)) { 

 1383  copy2>setInteger(iColumn); 

 1384  int kOrig = original[iColumn]; 

 1385  setOptionalInteger(kOrig); 

 1386  } 

 1387  } 

 1388  del[nDel++] = i; 

 1389  } 

 1390  } 

 1391  } 

 1392  } 

 1393  } 

 1394  if (nDel) { 

 1395  copy2>deleteRows(nDel, del); 

 1396  } 

 1397  if (nDel != numberRows) { 

 1398  nDel = 0; 

 1399  for (int i = 0; i < numberColumns; i++) { 

 1400  if (copy2>isInteger(i)) { 

 1401  del[nDel++] = i; 

 1402  nExtra++; 

 1403  } 

 1404  } 

 1405  } else { 

 1406  nDel = 0; 

 1407  } 

 1408  delete copy1; 

 1409  copy1 = copy2>clone(); 

 1410  delete copy2; 

 1411  } 

 1412  // See if what's left is a network 

 1413  bool couldBeNetwork = false; 

 1414  if (copy1>getNumRows() && copy1>getNumCols()) { 

 1415  #ifdef COIN_HAS_CLP 

 1416  OsiClpSolverInterface * clpSolver 

 1417  = dynamic_cast<OsiClpSolverInterface *> (copy1); 

 1418  if (false && clpSolver) { 

 1419  numberRows = clpSolver>getNumRows(); 

 1420  char * rotate = new char[numberRows]; 

 1421  int n = clpSolver>getModelPtr()>findNetwork(rotate, 1.0); 

 1422  delete [] rotate; 

 1423  #ifdef CLP_INVESTIGATE 

 1424  printf("INTA network %d rows out of %d\n", n, numberRows); 

 1425  #endif 

 1426  if (CoinAbs(n) == numberRows) { 

 1427  couldBeNetwork = true; 

 1428  for (int i = 0; i < numberRows; i++) { 

 1429  if (!possibleRow[i]) { 

 1430  couldBeNetwork = false; 

 1431  #ifdef CLP_INVESTIGATE 

 1432  printf("but row %d is bad\n", i); 

 1433  #endif 

 1434  break; 

 1435  } 

 1436  } 

 1437  } 

 1438  } else 

 1439  #endif 

 1440  { 

 1441  numberColumns = copy1>getNumCols(); 

 1442  numberRows = copy1>getNumRows(); 

 1443  const double * rowLower = copy1>getRowLower(); 

 1444  const double * rowUpper = copy1>getRowUpper(); 

 1445  couldBeNetwork = true; 

 1446  for (int i = 0; i < numberRows; i++) { 

 1447  if (rowLower[i] > 1.0e20 && fabs(rowLower[i]  floor(rowLower[i] + 0.5)) > 1.0e12) { 

 1448  couldBeNetwork = false; 

 1449  break; 

 1450  } 

 1451  if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i]  floor(rowUpper[i] + 0.5)) > 1.0e12) { 

 1452  couldBeNetwork = false; 

 1453  break; 

 1454  } 

 1455  } 

 1456  if (couldBeNetwork) { 

 1457  const CoinPackedMatrix * matrixByCol = copy1>getMatrixByCol(); 

 1458  const double * element = matrixByCol>getElements(); 

 1459  //const int * row = matrixByCol>getIndices(); 

 1460  const CoinBigIndex * columnStart = matrixByCol>getVectorStarts(); 

 1461  const int * columnLength = matrixByCol>getVectorLengths(); 

 1462  for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 

 1463  CoinBigIndex start = columnStart[iColumn]; 

 1464  CoinBigIndex end = start + columnLength[iColumn]; 

 1465  if (end > start + 2) { 

 1466  couldBeNetwork = false; 

 1467  break; 

 1468  } 

 1469  int type = 0; 

 1470  for (CoinBigIndex j = start; j < end; j++) { 

 1471  double value = element[j]; 

 1472  if (fabs(value) != 1.0) { 

 1473  couldBeNetwork = false; 

 1474  break; 

 1475  } else if (value == 1.0) { 

 1476  if ((type&1) == 0) 

 1477  type = 1; 

 1478  else 

 1479  type = 7; 

 1480  } else if (value == 1.0) { 

 1481  if ((type&2) == 0) 

 1482  type = 2; 

 1483  else 

 1484  type = 7; 

 1485  } 

 1486  } 

 1487  if (type > 3) { 

 1488  couldBeNetwork = false; 

 1489  break; 

 1490  } 

 1491  } 

 1492  } 

 1493  } 

 1494  } 

 1495  if (couldBeNetwork) { 

 1496  for (int i = 0; i < numberColumns; i++) 

 1497  setOptionalInteger(original[i]); 

 1498  } 

 1499  if (nExtra  couldBeNetwork) { 

 1500  numberColumns = copy1>getNumCols(); 

 1501  numberRows = copy1>getNumRows(); 

 1502  if (!numberColumns  !numberRows) { 

 1503  int numberColumns = solver_>getNumCols(); 

 1504  for (int i = 0; i < numberColumns; i++) 

 1505  assert(solver_>isInteger(i)); 

 1506  } 

 1507  #ifdef CLP_INVESTIGATE 

 1508  if (couldBeNetwork  nExtra) 

 1509  printf("INTA %d extra integers, %d left%s\n", nExtra, 

 1510  numberColumns, 

 1511  couldBeNetwork ? ", all network" : ""); 

 1512  #endif 

 1513  findIntegers(true, 2); 

 1514  convertToDynamic(); 

 1515  } 

 1516  #ifdef CLP_INVESTIGATE 

 1517  if (!couldBeNetwork && copy1>getNumCols() && 

 1518  copy1>getNumRows()) { 

 1519  printf("INTA %d rows and %d columns remain\n", 

 1520  copy1>getNumRows(), copy1>getNumCols()); 

 1521  if (copy1>getNumCols() < 200) { 

 1522  copy1>writeMps("moreint"); 

 1523  printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n", 

 1524  copy1>getNumRows(), copy1>getNumCols()); 

 1525  } 

 1526  } 

 1527  #endif 

 1528  delete copy1; 

 1529  delete [] del; 

 1530  delete [] original; 

 1531  delete [] possibleRow; 

 1532  // double check increment 

 1533  analyzeObjective(); 

 1534  } 

[2]  1535  /** 

 1536  \todo 

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

 1538  routine, after calling the solver for an initial solution 

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

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

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

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

 1543  checked. 

 1544  

 1545  */ 

 1546  

 1547  /* 

 1548  The overall flow can be divided into three stages: 

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

 1550  do all the setup for B&C. 

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

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

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

[1286]  1554  

[2]  1555  Keep in mind that in general there is no node in the search tree that 

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

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

 1558  system held by the solver. 

 1559  */ 

[1132]  1560  #ifdef COIN_HAS_CPX 

 1561  #include "OsiCpxSolverInterface.hpp" 

[1271]  1562  #include "cplex.h" 

[1132]  1563  #endif 

[1286]  1564  void CbcModel::branchAndBound(int doStatistics) 

[2]  1565  

 1566  { 

[1286]  1567  /* 

[1627]  1568  Capture a time stamp before we start (unless set). 

[1286]  1569  */ 

[1627]  1570  if (!dblParam_[CbcStartSeconds]) { 

 1571  if (!useElapsedTime()) 

 1572  dblParam_[CbcStartSeconds] = CoinCpuTime(); 

 1573  else 

 1574  dblParam_[CbcStartSeconds] = CoinGetTimeOfDay(); 

 1575  } 

[1286]  1576  dblParam_[CbcSmallestChange] = COIN_DBL_MAX; 

 1577  dblParam_[CbcSumChange] = 0.0; 

 1578  dblParam_[CbcLargestChange] = 0.0; 

 1579  intParam_[CbcNumberBranches] = 0; 

[1315]  1580  // Double check optimization directions line up 

 1581  dblParam_[CbcOptimizationDirection] = solver_>getObjSense(); 

[1286]  1582  strongInfo_[0] = 0; 

 1583  strongInfo_[1] = 0; 

 1584  strongInfo_[2] = 0; 

 1585  strongInfo_[3] = 0; 

 1586  strongInfo_[4] = 0; 

 1587  strongInfo_[5] = 0; 

 1588  strongInfo_[6] = 0; 

 1589  numberStrongIterations_ = 0; 

 1590  currentNode_ = NULL; 

 1591  // See if should do cuts old way 

 1592  if (parallelMode() < 0) { 

 1593  specialOptions_ = 4096 + 8192; 

 1594  } else if (parallelMode() > 0) { 

 1595  specialOptions_ = 4096; 

 1596  } 

 1597  int saveMoreSpecialOptions = moreSpecialOptions_; 

 1598  if (dynamic_cast<CbcTreeLocal *> (tree_)) 

 1599  specialOptions_ = 4096 + 8192; 

[838]  1600  #ifdef COIN_HAS_CLP 

[1286]  1601  { 

 1602  OsiClpSolverInterface * clpSolver 

 1603  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 1604  if (clpSolver) { 

 1605  // pass in disaster handler 

 1606  CbcDisasterHandler handler(this); 

 1607  clpSolver>passInDisasterHandler(&handler); 

 1608  // Initialise solvers seed 

 1609  clpSolver>getModelPtr()>setRandomSeed(1234567); 

[1393]  1610  #ifdef JJF_ZERO 

[1286]  1611  // reduce factorization frequency 

 1612  int frequency = clpSolver>getModelPtr()>factorizationFrequency(); 

 1613  clpSolver>getModelPtr()>setFactorizationFrequency(CoinMin(frequency, 120)); 

[1016]  1614  #endif 

[1286]  1615  } 

 1616  } 

[838]  1617  #endif 

[1286]  1618  // original solver (only set if preprocessing) 

 1619  OsiSolverInterface * originalSolver = NULL; 

 1620  int numberOriginalObjects = numberObjects_; 

 1621  OsiObject ** originalObject = NULL; 

 1622  // Save whether there were any objects 

 1623  bool noObjects = (numberObjects_ == 0); 

 1624  // Set up strategies 

[1409]  1625  /* 

 1626  See if the user has supplied a strategy object and deal with it if present. 

 1627  The call to setupOther will set numberStrong_ and numberBeforeTrust_, and 

 1628  perform integer preprocessing, if requested. 

[1364]  1629  

[1409]  1630  We need to hang on to a pointer to solver_. setupOther will assign a 

 1631  preprocessed solver to model, but will instruct assignSolver not to trash the 

 1632  existing one. 

 1633  */ 

[1357]  1634  if (strategy_) { 

 1635  // May do preprocessing 

 1636  originalSolver = solver_; 

 1637  strategy_>setupOther(*this); 

 1638  if (strategy_>preProcessState()) { 

 1639  // preprocessing done 

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

 1641  // infeasible 

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

 1643  status_ = 0 ; 

 1644  secondaryStatus_ = 1; 

 1645  originalContinuousObjective_ = COIN_DBL_MAX; 

 1646  return ; 

 1647  } else if (numberObjects_ && object_) { 

 1648  numberOriginalObjects = numberObjects_; 

 1649  // redo sequence 

 1650  numberIntegers_ = 0; 

 1651  int numberColumns = getNumCols(); 

 1652  int nOrig = originalSolver>getNumCols(); 

 1653  CglPreProcess * process = strategy_>process(); 

 1654  assert (process); 

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

 1656  // allow for cliques etc 

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

 1658  // try and redo debugger 

 1659  OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_>getRowCutDebuggerAlways()); 

 1660  if (debugger) 

 1661  debugger>redoSolution(numberColumns, originalColumns); 

[1409]  1662  // Userprovided solution might have been best. Synchronise. 

[1357]  1663  if (bestSolution_) { 

 1664  // need to redo  in case no better found in BAB 

 1665  // just get integer part right 

 1666  for (int i = 0; i < numberColumns; i++) { 

 1667  int jColumn = originalColumns[i]; 

 1668  bestSolution_[i] = bestSolution_[jColumn]; 

 1669  } 

 1670  } 

 1671  originalObject = object_; 

 1672  // object number or 1 

 1673  int * temp = new int[nOrig]; 

 1674  int iColumn; 

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

 1676  temp[iColumn] = 1; 

 1677  int iObject; 

 1678  int nNonInt = 0; 

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

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

 1681  if (iColumn < 0) { 

 1682  nNonInt++; 

 1683  } else { 

 1684  temp[iColumn] = iObject; 

 1685  } 

 1686  } 

 1687  int numberNewIntegers = 0; 

 1688  int numberOldIntegers = 0; 

 1689  int numberOldOther = 0; 

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

 1691  int jColumn = originalColumns[iColumn]; 

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

 1693  int iObject = temp[jColumn]; 

 1694  CbcSimpleInteger * obj = 

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

 1696  if (obj) 

 1697  numberOldIntegers++; 

 1698  else 

 1699  numberOldOther++; 

 1700  } else if (isInteger(iColumn)) { 

 1701  numberNewIntegers++; 

 1702  } 

 1703  } 

 1704  /* 

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

 1706  Make a large enough array for all objects 

 1707  */ 

 1708  numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt; 

 1709  object_ = new OsiObject * [numberObjects_]; 

 1710  delete [] integerVariable_; 

 1711  integerVariable_ = new int [numberNewIntegers+numberOldIntegers]; 

 1712  /* 

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

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

 1715  lower bounds. 

 1716  */ 

 1717  numberIntegers_ = 0; 

 1718  int n = originalColumns[numberColumns1] + 1; 

 1719  int * backward = new int[n]; 

 1720  int i; 

 1721  for ( i = 0; i < n; i++) 

 1722  backward[i] = 1; 

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

 1724  backward[originalColumns[i]] = i; 

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

 1726  int jColumn = originalColumns[iColumn]; 

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

 1728  int iObject = temp[jColumn]; 

 1729  CbcSimpleInteger * obj = 

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

 1731  if (obj) { 

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

 1733  // redo ids etc 

 1734  //object_[numberIntegers_]>resetSequenceEtc(numberColumns,originalColumns); 

 1735  object_[numberIntegers_]>resetSequenceEtc(numberColumns, backward); 

 1736  integerVariable_[numberIntegers_++] = iColumn; 

 1737  } 

 1738  } else if (isInteger(iColumn)) { 

 1739  object_[numberIntegers_] = 

 1740  new CbcSimpleInteger(this, iColumn); 

 1741  integerVariable_[numberIntegers_++] = iColumn; 

 1742  } 

 1743  } 

 1744  delete [] backward; 

 1745  numberObjects_ = numberIntegers_; 

 1746  // Now append other column stuff 

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

 1748  int jColumn = originalColumns[iColumn]; 

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

 1750  int iObject = temp[jColumn]; 

 1751  CbcSimpleInteger * obj = 

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

 1753  if (!obj) { 

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

 1755  // redo ids etc 

 1756  CbcObject * obj = 

 1757  dynamic_cast <CbcObject *>(object_[numberObjects_]) ; 

 1758  assert (obj); 

 1759  obj>redoSequenceEtc(this, numberColumns, originalColumns); 

 1760  numberObjects_++; 

 1761  } 

 1762  } 

 1763  } 

 1764  // now append non column stuff 

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

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

 1767  if (iColumn < 0) { 

 1768  // already has column numbers changed 

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

[1393]  1770  #ifdef JJF_ZERO 

[1357]  1771  // redo ids etc 

 1772  CbcObject * obj = 

 1773  dynamic_cast <CbcObject *>(object_[numberObjects_]) ; 

 1774  assert (obj); 

 1775  obj>redoSequenceEtc(this, numberColumns, originalColumns); 

 1776  #endif 

 1777  numberObjects_++; 

 1778  } 

 1779  } 

 1780  delete [] temp; 

 1781  if (!numberObjects_) 

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

 1783  } else { 

 1784  int numberColumns = getNumCols(); 

 1785  CglPreProcess * process = strategy_>process(); 

 1786  assert (process); 

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

 1788  // try and redo debugger 

 1789  OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_>getRowCutDebuggerAlways()); 

 1790  if (debugger) 

 1791  debugger>redoSolution(numberColumns, originalColumns); 

 1792  } 

 1793  } else { 

 1794  //no preprocessing 

 1795  originalSolver = NULL; 

 1796  } 

 1797  strategy_>setupCutGenerators(*this); 

 1798  strategy_>setupHeuristics(*this); 

 1799  // Set strategy print level to models 

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

 1801  } 

[1286]  1802  eventHappened_ = false; 

 1803  CbcEventHandler *eventHandler = getEventHandler() ; 

 1804  if (eventHandler) 

 1805  eventHandler>setModel(this); 

[1016]  1806  #define CLIQUE_ANALYSIS 

[833]  1807  #ifdef CLIQUE_ANALYSIS 

[1286]  1808  // set up for probing 

[1387]  1809  // If we're doing clever stuff with cliques, additional info here. 

[1286]  1810  if (!parentModel_) 

 1811  probingInfo_ = new CglTreeProbingInfo(solver_); 

 1812  else 

 1813  probingInfo_ = NULL; 

[833]  1814  #else 

[1286]  1815  probingInfo_ = NULL; 

[833]  1816  #endif 

[640]  1817  

[1286]  1818  // Try for dominated columns 

 1819  if ((specialOptions_&64) != 0) { 

 1820  CglDuplicateRow dupcuts(solver_); 

 1821  dupcuts.setMode(2); 

 1822  CglStored * storedCuts = dupcuts.outDuplicates(solver_); 

 1823  if (storedCuts) { 

 1824  printf("adding dup cuts\n"); 

 1825  addCutGenerator(storedCuts, 1, "StoredCuts from dominated", 

 1826  true, false, false, 200); 

 1827  } 

[1271]  1828  } 

[1286]  1829  if (!nodeCompare_) 

 1830  nodeCompare_ = new CbcCompareDefault();; 

 1831  // See if hot start wanted 

 1832  CbcCompareBase * saveCompare = NULL; 

[1387]  1833  // User supplied hotstart. Adapt for preprocessing. 

[1286]  1834  if (hotstartSolution_) { 

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

 1836  CglPreProcess * process = strategy_>process(); 

 1837  assert (process); 

 1838  int n = solver_>getNumCols(); 

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

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

 1841  double * tempS = new double[n]; 

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

 1843  int iColumn = originalColumns[i]; 

 1844  tempS[i] = hotstartSolution_[iColumn]; 

 1845  } 

 1846  delete [] hotstartSolution_; 

 1847  hotstartSolution_ = tempS; 

 1848  if (hotstartPriorities_) { 

 1849  int * tempP = new int [n]; 

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

 1851  int iColumn = originalColumns[i]; 

 1852  tempP[i] = hotstartPriorities_[iColumn]; 

 1853  } 

 1854  delete [] hotstartPriorities_; 

 1855  hotstartPriorities_ = tempP; 

 1856  } 

[231]  1857  } 

[1286]  1858  saveCompare = nodeCompare_; 

 1859  // depth first 

 1860  nodeCompare_ = new CbcCompareDepth(); 

[231]  1861  } 

[1286]  1862  if (!problemFeasibility_) 

 1863  problemFeasibility_ = new CbcFeasibilityBase(); 

[2]  1864  # ifdef CBC_DEBUG 

[1286]  1865  std::string problemName ; 

 1866  solver_>getStrParam(OsiProbName, problemName) ; 

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

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

[2]  1869  # endif 

[1286]  1870  /* 

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

 1872  */ 

 1873  status_ = 0 ; 

 1874  secondaryStatus_ = 0; 

 1875  phase_ = 0; 

 1876  /* 

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

 1878  CbcSimpleInteger object for each integer variable. 

 1879  */ 

 1880  findIntegers(false) ; 

 1881  // Say not dynamic pseudo costs 

 1882  ownership_ &= ~0x40000000; 

 1883  // If dynamic pseudo costs then do 

 1884  if (numberBeforeTrust_) 

 1885  convertToDynamic(); 

[1387]  1886  // Set up char array to say if integer (speed) 

[1286]  1887  delete [] integerInfo_; 

 1888  { 

 1889  int n = solver_>getNumCols(); 

 1890  integerInfo_ = new char [n]; 

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

 1892  if (solver_>isInteger(i)) 

 1893  integerInfo_[i] = 1; 

 1894  else 

 1895  integerInfo_[i] = 0; 

 1896  } 

[197]  1897  } 

[1286]  1898  if (preferredWay_) { 

 1899  // set all unset ones 

 1900  for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) { 

 1901  CbcObject * obj = 

 1902  dynamic_cast <CbcObject *>(object_[iObject]) ; 

 1903  if (obj && !obj>preferredWay()) 

 1904  obj>setPreferredWay(preferredWay_); 

 1905  } 

[640]  1906  } 

[1286]  1907  /* 

 1908  Ensure that objects on the lists of OsiObjects, heuristics, and cut 

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

 1910  */ 

 1911  synchronizeModel() ; 

 1912  if (!solverCharacteristics_) { 

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

 1914  if (solverCharacteristics) { 

 1915  solverCharacteristics_ = solverCharacteristics; 

 1916  } else { 

 1917  // replace in solver 

 1918  OsiBabSolver defaultC; 

 1919  solver_>setAuxiliaryInfo(&defaultC); 

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

 1921  } 

[395]  1922  } 

[862]  1923  

[1286]  1924  solverCharacteristics_>setSolver(solver_); 

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

 1926  continuousObjective_ = COIN_DBL_MAX ; 

 1927  /* 

 1928  Solve the relaxation. 

[2]  1929  

[1286]  1930  Apparently there are circumstances where this will be nontrivial  i.e., 

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

 1932  continuous relaxation. 

 1933  */ 

 1934  /* Tell solver we are in Branch and Cut 

 1935  Could use last parameter for subtle differences */ 

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

[1053]  1937  #ifdef COIN_HAS_CLP 

[1286]  1938  { 

 1939  OsiClpSolverInterface * clpSolver 

 1940  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 1941  if (clpSolver) { 

 1942  ClpSimplex * clpSimplex = clpSolver>getModelPtr(); 

 1943  if ((specialOptions_&32) == 0) { 

 1944  // take off names 

 1945  clpSimplex>dropNames(); 

 1946  } 

 1947  // no crunch if mostly continuous 

 1948  if ((clpSolver>specialOptions()&(1 + 8)) != (1 + 8)) { 

 1949  int numberColumns = solver_>getNumCols(); 

 1950  if (numberColumns > 1000 && numberIntegers_*4 < numberColumns) 

 1951  clpSolver>setSpecialOptions(clpSolver>specialOptions()&(~1)); 

 1952  } 

 1953  //#define NO_CRUNCH 

[1053]  1954  #ifdef NO_CRUNCH 

[1286]  1955  printf("TEMP switching off crunch\n"); 

 1956  int iOpt = clpSolver>specialOptions(); 

 1957  iOpt &= ~1; 

 1958  iOpt = 65536; 

 1959  clpSolver>setSpecialOptions(iOpt); 

[1053]  1960  #endif 

[1286]  1961  } 

[1053]  1962  } 

 1963  #endif 

[1286]  1964  bool feasible; 

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

 1966  if (!solverCharacteristics_>solverType()  solverCharacteristics_>solverType() == 4) { 

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

[480]  1968  } else { 

[1286]  1969  // pick up given status 

 1970  feasible = (solver_>isProvenOptimal() && 

 1971  !solver_>isDualObjectiveLimitReached()) ; 

[480]  1972  } 

[1286]  1973  if (problemFeasibility_>feasible(this, 0) < 0) { 

 1974  feasible = false; // pretend infeasible 

[640]  1975  } 

[1286]  1976  numberSavedSolutions_ = 0; 

 1977  int saveNumberStrong = numberStrong_; 

 1978  int saveNumberBeforeTrust = numberBeforeTrust_; 

 1979  /* 

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

 1981  continue with processing the root node. 

 1982  */ 

 1983  if (!feasible) { 

 1984  status_ = 0 ; 

 1985  if (!solver_>isProvenDualInfeasible()) { 

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

 1987  secondaryStatus_ = 1; 

 1988  } else { 

 1989  handler_>message(CBC_UNBOUNDED, messages_) << CoinMessageEol ; 

 1990  secondaryStatus_ = 7; 

 1991  } 

 1992  originalContinuousObjective_ = COIN_DBL_MAX; 

 1993  solverCharacteristics_ = NULL; 

 1994  return ; 

 1995  } else if (!numberObjects_) { 

 1996  // nothing to do 

 1997  solverCharacteristics_ = NULL; 

 1998  bestObjective_ = solver_>getObjValue() * solver_>getObjSense(); 

 1999  int numberColumns = solver_>getNumCols(); 

 2000  delete [] bestSolution_; 

 2001  bestSolution_ = new double[numberColumns]; 

 2002  CoinCopyN(solver_>getColSolution(), numberColumns, bestSolution_); 

 2003  return ; 

[640]  2004  } 

[1409]  2005  /* 

 2006  See if we're using the Osi side of the branching hierarchy. If so, either 

 2007  convert existing CbcObjects to OsiObjects, or generate them fresh. In the 

 2008  first case, CbcModel owns the objects on the object_ list. In the second 

 2009  case, the solver holds the objects and object_ simply points to the 

 2010  solver's list. 

[1364]  2011  

[1409]  2012  080417 The conversion code here (the block protected by `if (obj)') cannot 

 2013  possibly be correct. On the Osi side, descent is OsiObject > OsiObject2 > 

 2014  all other Osi object classes. On the Cbc side, it's OsiObject > CbcObject 

 2015  > all other Cbc object classes. It's structurally impossible for any Osi 

 2016  object to descend from CbcObject. The only thing I can see is that this is 

 2017  really dead code, and object detection is now handled from the Osi side. 

 2018  */ 

[1286]  2019  // Convert to Osi if wanted 

 2020  bool useOsiBranching = false; 

 2021  //OsiBranchingInformation * persistentInfo = NULL; 

 2022  if (branchingMethod_ && branchingMethod_>chooseMethod()) { 

 2023  useOsiBranching = true; 

 2024  //persistentInfo = new OsiBranchingInformation(solver_); 

 2025  if (numberOriginalObjects) { 

 2026  for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) { 

 2027  CbcObject * obj = 

 2028  dynamic_cast <CbcObject *>(object_[iObject]) ; 

 2029  if (obj) { 

 2030  CbcSimpleInteger * obj2 = 

 2031  dynamic_cast <CbcSimpleInteger *>(obj) ; 

 2032  if (obj2) { 

 2033  // back to Osi land 

 2034  object_[iObject] = obj2>osiObject(); 

 2035  delete obj; 

 2036  } else { 

 2037  OsiSimpleInteger * obj3 = 

 2038  dynamic_cast <OsiSimpleInteger *>(obj) ; 

 2039  if (!obj3) { 

 2040  OsiSOS * obj4 = 

 2041  dynamic_cast <OsiSOS *>(obj) ; 

 2042  if (!obj4) { 

 2043  CbcSOS * obj5 = 

 2044  dynamic_cast <CbcSOS *>(obj) ; 

 2045  if (obj5) { 

 2046  // back to Osi land 

 2047  object_[iObject] = obj5>osiObject(solver_); 

 2048  } else { 

 2049  printf("Code up CbcObject type in Osi land\n"); 

 2050  abort(); 

 2051  } 

 2052  } 

 2053  } 

 2054  } 

 2055  } 

 2056  } 

 2057  // and add to solver 

 2058  //if (!solver_>numberObjects()) { 

 2059  solver_>addObjects(numberObjects_, object_); 

 2060  //} else { 

 2061  //if (solver_>numberObjects()!=numberOriginalObjects) { 

 2062  //printf("should have trapped that solver has objects before\n"); 

 2063  //abort(); 

 2064  //} 

 2065  //} 

 2066  } else { 

[1409]  2067  /* 

 2068  As of 080104, findIntegersAndSOS is misleading  the default OSI 

 2069  implementation finds only integers. 

 2070  */ 

[1286]  2071  // do from solver 

 2072  deleteObjects(false); 

 2073  solver_>findIntegersAndSOS(false); 

 2074  numberObjects_ = solver_>numberObjects(); 

 2075  object_ = solver_>objects(); 

 2076  ownObjects_ = false; 

 2077  } 

 2078  branchingMethod_>chooseMethod()>setSolver(solver_); 

[640]  2079  } 

[1387]  2080  // take off heuristics if have to (some do not work with SOS, for example) 

 2081  // object should know what's safe. 

[1286]  2082  { 

 2083  int numberOdd = 0; 

 2084  int numberSOS = 0; 

 2085  for (int i = 0; i < numberObjects_; i++) { 

 2086  if (!object_[i]>canDoHeuristics()) 

 2087  numberOdd++; 

 2088  CbcSOS * obj = 

 2089  dynamic_cast <CbcSOS *>(object_[i]) ; 

 2090  if (obj) 

 2091  numberSOS++; 

 2092  } 

 2093  if (numberOdd) { 

 2094  if (numberHeuristics_) { 

 2095  int k = 0; 

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

 2097  if (!heuristic_[i]>canDealWithOdd()) 

 2098  delete heuristic_[i]; 

 2099  else 

 2100  heuristic_[k++] = heuristic_[i]; 

 2101  } 

 2102  if (!k) { 

 2103  delete [] heuristic_; 

 2104  heuristic_ = NULL; 

 2105  } 

 2106  numberHeuristics_ = k; 

 2107  handler_>message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ; 

 2108  } 

 2109  } else if (numberSOS) { 

 2110  specialOptions_ = 128; // say can do SOS in dynamic mode 

 2111  // switch off fast nodes for now 

 2112  fastNodeDepth_ = 1; 

 2113  } 

 2114  if (numberThreads_ > 0) { 

 2115  // switch off fast nodes for now 

 2116  fastNodeDepth_ = 1; 

 2117  } 

[931]  2118  } 

[1286]  2119  // Save objective (just so user can access it) 

 2120  originalContinuousObjective_ = solver_>getObjValue(); 

 2121  bestPossibleObjective_ = originalContinuousObjective_; 

 2122  sumChangeObjective1_ = 0.0; 

 2123  sumChangeObjective2_ = 0.0; 

 2124  /* 

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

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

 2127  see if it cuts off the optimal answer. 

 2128  */ 

 2129  // If debugger exists set specialOptions_ bit 

 2130  if (solver_>getRowCutDebuggerAlways()) { 

 2131  specialOptions_ = 1; 

 2132  } 

[56]  2133  

 2134  # ifdef CBC_DEBUG 

[1286]  2135  if ((specialOptions_&1) == 0) 

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

 2137  if (solver_>getRowCutDebuggerAlways()) 

 2138  specialOptions_ = 1; 

[2]  2139  # endif 

 2140  

[1286]  2141  /* 

 2142  Begin setup to process a feasible root node. 

 2143  */ 

 2144  bestObjective_ = CoinMin(bestObjective_, 1.0e50) ; 

 2145  if (!bestSolution_) { 

 2146  numberSolutions_ = 0 ; 

 2147  numberHeuristicSolutions_ = 0 ; 

 2148  } 

 2149  stateOfSearch_ = 0; 

 2150  // Everything is minimization 

 2151  { 

 2152  // needed to sync cutoffs 

 2153  double value ; 

 2154  solver_>getDblParam(OsiDualObjectiveLimit, value) ; 

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

 2156  } 

 2157  double cutoff = getCutoff() ; 

 2158  double direction = solver_>getObjSense() ; 

 2159  dblParam_[CbcOptimizationDirection] = direction; 

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

 2161  messageHandler()>message(CBC_CUTOFF_WARNING1, 

 2162  messages()) 

 2163  << cutoff << cutoff << CoinMessageEol ; 

 2164  if (cutoff > bestObjective_) 

 2165  cutoff = bestObjective_ ; 

 2166  setCutoff(cutoff) ; 

 2167  /* 

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

 2169  */ 

 2170  int numberColumns = getNumCols() ; 

 2171  if (!currentSolution_) 

 2172  currentSolution_ = new double[numberColumns] ; 

 2173  testSolution_ = currentSolution_; 

 2174  /* 

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

 2176  constraint system (aka the continuous system). 

 2177  */ 

 2178  continuousSolver_ = solver_>clone() ; 

[152]  2179  

[1286]  2180  numberRowsAtContinuous_ = getNumRows() ; 

 2181  solver_>saveBaseModel(); 

 2182  /* 

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

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

 2185  installed. 

 2186  */ 

 2187  if (solverCharacteristics_>reducedCostsAccurate()) 

 2188  analyzeObjective() ; 

 2189  { 

 2190  // may be able to change cutoff now 

 2191  double cutoff = getCutoff(); 

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

 2193  if (cutoff > bestObjective_  increment) { 

 2194  cutoff = bestObjective_  increment ; 

 2195  setCutoff(cutoff) ; 

 2196  } 

[1271]  2197  } 

[931]  2198  #ifdef COIN_HAS_CLP 

[1286]  2199  // Possible save of pivot method 

 2200  ClpDualRowPivot * savePivotMethod = NULL; 

 2201  { 

 2202  // pass tolerance and increment to solver 

 2203  OsiClpSolverInterface * clpSolver 

 2204  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 2205  if (clpSolver) 

 2206  clpSolver>setStuff(getIntegerTolerance(), getCutoffIncrement()); 

 2207  } 

[931]  2208  #endif 

[1286]  2209  /* 

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

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

 2212  that produced a given cut. 

 2213  */ 

 2214  maximumWhich_ = 1000 ; 

 2215  delete [] whichGenerator_; 

 2216  whichGenerator_ = new int[maximumWhich_] ; 

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

 2218  maximumNumberCuts_ = 0 ; 

 2219  currentNumberCuts_ = 0 ; 

 2220  delete [] addedCuts_ ; 

 2221  addedCuts_ = NULL ; 

 2222  OsiObject ** saveObjects = NULL; 

 2223  maximumRows_ = numberRowsAtContinuous_; 

 2224  currentDepth_ = 0; 

 2225  workingBasis_.resize(maximumRows_, numberColumns); 

 2226  /* 

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

 2228  (problems which require further exploration). 

 2229  */ 

 2230  CbcCompareDefault * compareActual 

 2231  = dynamic_cast<CbcCompareDefault *> (nodeCompare_); 

 2232  if (compareActual) { 

 2233  compareActual>setBestPossible(direction*solver_>getObjValue()); 

 2234  compareActual>setCutoff(getCutoff()); 

[1393]  2235  #ifdef JJF_ZERO 

[1357]  2236  if (false && !numberThreads_ && !parentModel_) { 

[1286]  2237  printf("CbcTreeArray ? threads ? parentArray\n"); 

 2238  // Setup new style tree 

 2239  delete tree_; 

 2240  tree_ = new CbcTreeArray(); 

 2241  } 

[1341]  2242  #endif 

[1132]  2243  } 

[1286]  2244  tree_>setComparison(*nodeCompare_) ; 

 2245  /* 

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

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

 2248  */ 

 2249  maximumDepth_ = 10 ; 

 2250  delete [] walkback_ ; 

 2251  walkback_ = new CbcNodeInfo * [maximumDepth_] ; 

 2252  lastDepth_ = 0; 

 2253  delete [] lastNodeInfo_ ; 

 2254  lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ; 

 2255  delete [] lastNumberCuts_ ; 

 2256  lastNumberCuts_ = new int [maximumDepth_] ; 

 2257  maximumCuts_ = 100; 

 2258  lastNumberCuts2_ = 0; 

 2259  delete [] lastCut_; 

 2260  lastCut_ = new const OsiRowCut * [maximumCuts_]; 

 2261  /* 

 2262  Used to generate bound edits for CbcPartialNodeInfo. 

 2263  */ 

 2264  double * lowerBefore = new double [numberColumns] ; 

 2265  double * upperBefore = new double [numberColumns] ; 

 2266  /* 

[1409]  2267  Set up to run heuristics and generate cuts at the root node. The heavy 

 2268  lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts. 

[2]  2269  

[1409]  2270  To start, tell cut generators they can be a bit more aggressive at the 

 2271  root node. 

[1364]  2272  

[1409]  2273  QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve 

 2274  with cuts at root'. Is phase_ = 1 the correct indication when 

 2275  doHeurisiticsAtRoot is called to run heuristics outside of the main 

 2276  cut / heurisitc / reoptimise loop in solveWithCuts? 

[1364]  2277  

[1286]  2278  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy 

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

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

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

[2]  2282  

[1286]  2283  At the end of all this, any cuts will be recorded in cuts and also 

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

 2285  removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been 

 2286  adjusted accordingly). 

[2]  2287  

[1286]  2288  Tell cut generators they can be a bit more aggressive at root node 

 2289  

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

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

 2292  */ 

 2293  phase_ = 1; 

 2294  int iCutGenerator; 

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

[1409]  2296  // If parallel switch off global cuts 

 2297  if (numberThreads_) { 

 2298  generator_[iCutGenerator]>setGlobalCuts(false); 

 2299  generator_[iCutGenerator]>setGlobalCutsAtRoot(false); 

 2300  } 

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

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

[1271]  2303  } 

[1286]  2304  OsiCuts cuts ; 

 2305  int anyAction = 1 ; 

 2306  numberOldActiveCuts_ = 0 ; 

 2307  numberNewCuts_ = 0 ; 

 2308  // Array to mark solution 

 2309  delete [] usedInSolution_; 

 2310  usedInSolution_ = new int[numberColumns]; 

 2311  CoinZeroN(usedInSolution_, numberColumns); 

 2312  /* 

 2313  For printing totals and for CbcNode (numberNodes_) 

 2314  */ 

 2315  numberIterations_ = 0 ; 

 2316  numberSolves_ = 0 ; 

 2317  numberNodes_ = 0 ; 

 2318  numberNodes2_ = 0 ; 

 2319  maximumStatistics_ = 0; 

 2320  maximumDepthActual_ = 0; 

 2321  numberDJFixed_ = 0.0; 

[1315]  2322  if (!parentModel_) { 

 2323  if ((specialOptions_&262144) != 0) { 

 2324  // create empty stored cuts 

 2325  //storedRowCuts_ = new CglStored(solver_>getNumCols()); 

 2326  } else if ((specialOptions_&524288) != 0 && storedRowCuts_) { 

 2327  // tighten and set best solution 

 2328  // A) tight bounds on integer variables 

[1409]  2329  /* 

 2330  storedRowCuts_ are coming in from outside, probably for nonlinear. 

 2331  John was unsure about origin. 

 2332  */ 

[1315]  2333  const double * lower = solver_>getColLower(); 

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

 2335  const double * tightLower = storedRowCuts_>tightLower(); 

 2336  const double * tightUpper = storedRowCuts_>tightUpper(); 

 2337  int nTightened = 0; 

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

 2339  int iColumn = integerVariable_[i]; 

 2340  if (tightLower[iColumn] > lower[iColumn]) { 

 2341  nTightened++; 

 2342  solver_>setColLower(iColumn, tightLower[iColumn]); 

 2343  } 

 2344  if (tightUpper[iColumn] < upper[iColumn]) { 

 2345  nTightened++; 

 2346  solver_>setColUpper(iColumn, tightUpper[iColumn]); 

 2347  } 

 2348  } 

 2349  if (nTightened) 

 2350  printf("%d tightened by alternate cuts\n", nTightened); 

 2351  if (storedRowCuts_>bestObjective() < bestObjective_) { 

 2352  // B) best solution 

 2353  double objValue = storedRowCuts_>bestObjective(); 

 2354  setBestSolution(CBC_SOLUTION, objValue, 

 2355  storedRowCuts_>bestSolution()) ; 

 2356  // Do heuristics 

 2357  // Allow RINS 

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

 2359  CbcHeuristicRINS * rins 

 2360  = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]); 

 2361  if (rins) { 

 2362  rins>setLastNode(100); 

 2363  } 

 2364  } 

 2365  } 

 2366  } 

 2367  } 

[1409]  2368  /* 

 2369  Run heuristics at the root. This is the only opportunity to run FPump; it 

 2370  will be removed from the heuristics list by doHeuristicsAtRoot. 

 2371  */ 

[1286]  2372  // Do heuristics 

 2373  doHeuristicsAtRoot(); 

[1409]  2374  /* 

 2375  Grepping through the code, it would appear that this is a command line 

 2376  debugging hook. There's no obvious place in the code where this is set to 

 2377  a negative value. 

[1387]  2378  

[1409]  2379  User hook, says John. 

 2380  */ 

[1286]  2381  if ( intParam_[CbcMaxNumNode] < 0) 

 2382  eventHappened_ = true; // stop as fast as possible 

 2383  stoppedOnGap_ = false ; 

 2384  // See if can stop on gap 

 2385  bestPossibleObjective_ = solver_>getObjValue() * solver_>getObjSense(); 

 2386  double testGap = CoinMax(dblParam_[CbcAllowableGap], 

 2387  CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_)) 

 2388  * dblParam_[CbcAllowableFractionGap]); 

 2389  if (bestObjective_  bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) { 

 2390  if (bestPossibleObjective_ < getCutoff()) 

 2391  stoppedOnGap_ = true ; 

 2392  feasible = false; 

 2393  //eventHappened_=true; // stop as fast as possible 

[1271]  2394  } 

[1409]  2395  /* 

 2396  Set up for statistics collection, if requested. Standard values are 

 2397  documented in CbcModel.hpp. The magic number 100 will trigger a dump of 

 2398  CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another 

 2399  command line debugging hook. 

 2400  */ 

[1286]  2401  statistics_ = NULL; 

 2402  // Do on switch 

 2403  if (doStatistics > 0 && doStatistics <= 100) { 

 2404  maximumStatistics_ = 10000; 

 2405  statistics_ = new CbcStatistics * [maximumStatistics_]; 

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

[1271]  2407  } 

[1286]  2408  // See if we can add integers 

[1330]  2409  if (noObjects && numberIntegers_ < solver_>getNumCols() && (specialOptions_&65536) != 0 && !parentModel_) 

[1357]  2410  AddIntegers(); 

[1409]  2411  /* 

 2412  Do an initial round of cut generation for the root node. Depending on the 

 2413  type of underlying solver, we may want to do this even if the initial query 

 2414  to the objects indicates they're satisfied. 

[1330]  2415  

[1409]  2416  solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise 

 2417  loop (including reduced cost fixing) until no cuts are generated, the 

 2418  change in objective falls off, or the limit on the number of rounds of cut 

 2419  generation is exceeded. 

[1364]  2420  

[1409]  2421  At the end of all this, any cuts will be recorded in cuts and also 

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

 2423  removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been 

 2424  adjusted accordingly). 

 2425  */ 

[1330]  2426  int iObject ; 

 2427  int preferredWay ; 

 2428  int numberUnsatisfied = 0 ; 

 2429  delete [] currentSolution_; 

 2430  currentSolution_ = new double [numberColumns]; 

 2431  testSolution_ = currentSolution_; 

 2432  memcpy(currentSolution_, solver_>getColSolution(), 

 2433  numberColumns*sizeof(double)) ; 

 2434  // point to useful information 

 2435  OsiBranchingInformation usefulInfo = usefulInformation(); 

 2436  

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

 2438  double infeasibility = 

 2439  object_[iObject]>infeasibility(&usefulInfo, preferredWay) ; 

 2440  if (infeasibility ) numberUnsatisfied++ ; 

 2441  } 

 2442  // replace solverType 

 2443  if (solverCharacteristics_>tryCuts()) { 

 2444  

 2445  if (numberUnsatisfied) { 

 2446  // User event 

 2447  if (!eventHappened_ && feasible) { 

 2448  feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_, 

 2449  NULL); 

 2450  if ((specialOptions_&524288) != 0 && !parentModel_ 

 2451  && storedRowCuts_) { 

 2452  if (feasible) { 

 2453  /* pick up stuff and try again 

 2454  add cuts, maybe keep around 

 2455  do best solution and if so new heuristics 

 2456  obviously tighten bounds 

 2457  */ 

 2458  // A and B probably done on entry 

 2459  // A) tight bounds on integer variables 

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

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

 2462  const double * tightLower = storedRowCuts_>tightLower(); 

 2463  const double * tightUpper = storedRowCuts_>tightUpper(); 

 2464  int nTightened = 0; 

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

 2466  int iColumn = integerVariable_[i]; 

 2467  if (tightLower[iColumn] > lower[iColumn]) { 

 2468  nTightened++; 

 2469  solver_>setColLower(iColumn, tightLower[iColumn]); 

[1286]  2470  } 

[1330]  2471  if (tightUpper[iColumn] < upper[iColumn]) { 

 2472  nTightened++; 

 2473  solver_>setColUpper(iColumn, tightUpper[iColumn]); 

 2474  } 

[1286]  2475  } 

[1330]  2476  if (nTightened) 

 2477  printf("%d tightened by alternate cuts\n", nTightened); 

 2478  if (storedRowCuts_>bestObjective() < bestObjective_) { 

 2479  // B) best solution 

 2480  double objValue = storedRowCuts_>bestObjective(); 

 2481  setBestSolution(CBC_SOLUTION, objValue, 

 2482  storedRowCuts_>bestSolution()) ; 

 2483  // Do heuristics 

 2484  // Allow RINS 

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

 2486  CbcHeuristicRINS * rins 

 2487  = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]); 

 2488  if (rins) { 

 2489  rins>setLastNode(100); 

[1286]  2490  } 

 2491  } 

[1330]  2492  doHeuristicsAtRoot(); 

[1286]  2493  } 

[1393]  2494  #ifdef JJF_ZERO 

[1330]  2495  int nCuts = storedRowCuts_>sizeRowCuts(); 

 2496  // add to global list 

 2497  for (int i = 0; i < nCuts; i++) { 

 2498  OsiRowCut newCut(*storedRowCuts_>rowCutPointer(i)); 

 2499  newCut.setGloballyValidAsInteger(2); 

 2500  newCut.mutableRow().setTestForDuplicateIndex(false); 

 2501  globalCuts_.insert(newCut) ; 

[1286]  2502  } 

[1330]  2503  #else 

 2504  addCutGenerator(storedRowCuts_, 99, "Stored from previous run", 

 2505  true, false, false, 200); 

[1271]  2506  #endif 

[1330]  2507  // Set cuts as active 

 2508  delete [] addedCuts_ ; 

 2509  maximumNumberCuts_ = cuts.sizeRowCuts(); 

 2510  if (maximumNumberCuts_) { 

 2511  addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_]; 

 2512  } else { 

 2513  addedCuts_ = NULL; 

[1286]  2514  } 

[1330]  2515  for (int i = 0; i < maximumNumberCuts_; i++) 

 2516  addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i), 

 2517  NULL, 1, 1, 2); 

 2518  printf("size %d\n", cuts.sizeRowCuts()); 

 2519  cuts = OsiCuts(); 

 2520  currentNumberCuts_ = maximumNumberCuts_; 

 2521  feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_, 

 2522  NULL); 

 2523  for (int i = 0; i < maximumNumberCuts_; i++) 

 2524  delete addedCuts_[i]; 

[1286]  2525  } 

[1330]  2526  delete storedRowCuts_; 

 2527  storedRowCuts_ = NULL; 

[1286]  2528  } 

[1330]  2529  } else { 

 2530  feasible = false; 

[1286]  2531  } 

[1330]  2532  } else if (solverCharacteristics_>solutionAddsCuts()  

 2533  solverCharacteristics_>alwaysTryCutsAtRootNode()) { 

 2534  // may generate cuts and turn the solution 

 2535  //to an infeasible one 

 2536  feasible = solveWithCuts(cuts, 1, 

 2537  NULL); 

[1286]  2538  } 

[1271]  2539  } 

[1330]  2540  // check extra info on feasibility 

 2541  if (!solverCharacteristics_>mipFeasible()) 

 2542  feasible = false; 

[2]  2543  

[1286]  2544  // make cut generators less aggressive 

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

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

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

[1132]  2548  } 

[1286]  2549  currentNumberCuts_ = numberNewCuts_ ; 

 2550  // See if can stop on gap 

 2551  bestPossibleObjective_ = solver_>getObjValue() * solver_>getObjSense(); 

 2552  testGap = CoinMax(dblParam_[CbcAllowableGap], 

 2553  CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_)) 

 2554  * dblParam_[CbcAllowableFractionGap]); 

 2555  if (bestObjective_  bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) { 

 2556  if (bestPossibleObjective_ < getCutoff()) 

 2557  stoppedOnGap_ = true ; 

 2558  feasible = false; 

[1132]  2559  } 

[1286]  2560  // User event 

 2561  if (eventHappened_) 

 2562  feasible = false; 

 2563  #if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX) 

[1409]  2564  /* 

 2565  This is the notion of using Cbc stuff to get going, then calling cplex to 

 2566  finish off. 

 2567  */ 

[1286]  2568  if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == 2 && !parentModel_) { 

 2569  // Use Cplex to do search! 

 2570  double time1 = CoinCpuTime(); 

 2571  OsiClpSolverInterface * clpSolver 

 2572  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 2573  OsiCpxSolverInterface cpxSolver; 

 2574  double direction = clpSolver>getObjSense(); 

 2575  cpxSolver.setObjSense(direction); 

 2576  // load up cplex 

 2577  const CoinPackedMatrix * matrix = continuousSolver_>getMatrixByCol(); 

 2578  const double * rowLower = continuousSolver_>getRowLower(); 

 2579  const double * rowUpper = continuousSolver_>getRowUpper(); 

 2580  const double * columnLower = continuousSolver_>getColLower(); 

 2581  const double * columnUpper = continuousSolver_>getColUpper(); 

 2582  const double * objective = continuousSolver_>getObjCoefficients(); 

 2583  cpxSolver.loadProblem(*matrix, columnLower, columnUpper, 

 2584  objective, rowLower, rowUpper); 

 2585  double * setSol = new double [numberIntegers_]; 

 2586  int * setVar = new int [numberIntegers_]; 

 2587  // cplex doesn't know about objective offset 

 2588  double offset = clpSolver>getModelPtr()>objectiveOffset(); 

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

 2590  int iColumn = integerVariable_[i]; 

 2591  cpxSolver.setInteger(iColumn); 

 2592  if (bestSolution_) { 

 2593  setSol[i] = bestSolution_[iColumn]; 

 2594  setVar[i] = iColumn; 

 2595  } 

 2596  } 

 2597  CPXENVptr env = cpxSolver.getEnvironmentPtr(); 

 2598  CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL); 

 2599  cpxSolver.switchToMIP(); 

 2600  if (bestSolution_) { 

 2601  CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol); 

 2602  } 

 2603  if (clpSolver>getNumRows() > continuousSolver_>getNumRows() && false) { 

 2604  // add cuts 

 2605  const CoinPackedMatrix * matrix = clpSolver>getMatrixByRow(); 

 2606  const double * rhs = clpSolver>getRightHandSide(); 

 2607  const char * rowSense = clpSolver>getRowSense(); 

 2608  const double * elementByRow = matrix>getElements(); 

 2609  const int * column = matrix>getIndices(); 

 2610  const CoinBigIndex * rowStart = matrix>getVectorStarts(); 

 2611  const int * rowLength = matrix>getVectorLengths(); 

 2612  int nStart = continuousSolver_>getNumRows(); 

 2613  int nRows = clpSolver>getNumRows(); 

 2614  int size = rowStart[nRows1] + rowLength[nRows1]  

 2615  rowStart[nStart]; 

 2616  int nAdd = 0; 

 2617  double * rmatval = new double [size]; 

 2618  int * rmatind = new int [size]; 

 2619  int * rmatbeg = new int [nRowsnStart+1]; 

 2620  size = 0; 

 2621  rmatbeg[0] = 0; 

 2622  for (int i = nStart; i < nRows; i++) { 

 2623  for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) { 

 2624  rmatind[size] = column[k]; 

 2625  rmatval[size++] = elementByRow[k]; 

 2626  } 

 2627  nAdd++; 

 2628  rmatbeg[nAdd] = size; 

 2629  } 

 2630  CPXaddlazyconstraints(env, lpPtr, nAdd, size, 

 2631  rhs, rowSense, rmatbeg, 

 2632  rmatind, rmatval, NULL); 

 2633  CPXsetintparam( env, CPX_PARAM_REDUCE, 

 2634  // CPX_PREREDUCE_NOPRIMALORDUAL (0) 

 2635  CPX_PREREDUCE_PRIMALONLY); 

 2636  } 

 2637  if (getCutoff() < 1.0e50) { 

 2638  double useCutoff = getCutoff() + offset; 

 2639  if (bestObjective_ < 1.0e50) 

 2640  useCutoff = bestObjective_ + offset + 1.0e7; 

 2641  cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff* 

 2642  direction); 

 2643  if ( direction > 0.0 ) 

 2644  CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min 

 2645  else 

 2646  CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max 

 2647  } 

 2648  CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]); 

 2649  delete [] setSol; 

 2650  delete [] setVar; 

 2651  char printBuffer[200]; 

 2652  if (offset) { 

 2653  sprintf(printBuffer, "Add %g to all Cplex messages for true objective", 

 2654  offset); 

 2655  messageHandler()>message(CBC_GENERAL, messages()) 

 2656  << printBuffer << CoinMessageEol ; 

 2657  cpxSolver.setDblParam(OsiObjOffset, offset); 

 2658  } 

 2659  cpxSolver.branchAndBound(); 

 2660  double timeTaken = CoinCpuTime()  time1; 

 2661  sprintf(printBuffer, "Cplex took %g seconds", 

 2662  timeTaken); 

 2663  messageHandler()>message(CBC_GENERAL, messages()) 

 2664  << printBuffer << CoinMessageEol ; 

 2665  numberExtraNodes_ = CPXgetnodecnt(env, lpPtr); 

 2666  numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr); 

 2667  double value = cpxSolver.getObjValue() * direction; 

 2668  if (cpxSolver.isProvenOptimal() && value <= getCutoff()) { 

 2669  feasible = true; 

 2670  clpSolver>setWarmStart(NULL); 

 2671  // try and do solution 

 2672  double * newSolution = 

 2673  CoinCopyOfArray(cpxSolver.getColSolution(), 

 2674  getNumCols()); 

 2675  setBestSolution(CBC_STRONGSOL, value, newSolution) ; 

 2676  delete [] newSolution; 

 2677  } 

 2678  feasible = false; 

[1132]  2679  } 

[1286]  2680  #endif 

[1409]  2681  /* 

 2682  A hook to use clp to quickly explore some part of the tree. 

 2683  */ 

[1286]  2684  if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) { 

 2685  fastNodeDepth_ = 1; 

 2686  CbcObject * obj = 

 2687  new CbcFollowOn(this); 

 2688  obj>setPriority(1); 

 2689  addObjects(1, &obj); 

[1132]  2690  } 

[1286]  2691  int saveNumberSolves = numberSolves_; 

 2692  int saveNumberIterations = numberIterations_; 

 2693  if (fastNodeDepth_ >= 0 &&/*!parentModel_*/(specialOptions_&2048) == 0) { 

 2694  // add in a general depth object doClp 

 2695  int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : (fastNodeDepth_  100); 

 2696  CbcObject * obj = 

 2697  new CbcGeneralDepth(this, type); 

 2698  addObjects(1, &obj); 

 2699  // mark as done 

 2700  fastNodeDepth_ += 1000000; 

 2701  delete obj; 

 2702  // fake number of objects 

 2703  numberObjects_; 

 2704  if (parallelMode() < 1) { 

 2705  // But make sure position is correct 

 2706  OsiObject * obj2 = object_[numberObjects_]; 

 2707  obj = dynamic_cast<CbcObject *> (obj2); 

 2708  assert (obj); 

 2709  obj>setPosition(numberObjects_); 

 2710  } 

[1132]  2711  } 

[1286]  2712  #ifdef COIN_HAS_CLP 

 2713  #ifdef NO_CRUNCH 

 2714  if (true) { 

 2715  OsiClpSolverInterface * clpSolver 

 2716  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 2717  if (clpSolver && !parentModel_) { 

 2718  ClpSimplex * clpSimplex = clpSolver>getModelPtr(); 

 2719  clpSimplex>setSpecialOptions(clpSimplex>specialOptions()  131072); 

 2720  //clpSimplex>startPermanentArrays(); 

 2721  clpSimplex>setPersistenceFlag(2); 

 2722  } 

[1121]  2723  } 

[931]  2724  #endif 

 2725  #endif 

[1286]  2726  // Save copy of solver 

 2727  OsiSolverInterface * saveSolver = NULL; 

 2728  if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0) 

 2729  saveSolver = solver_>clone(); 

 2730  double checkCutoffForRestart = 1.0e100; 

[1330]  2731  saveModel(saveSolver, &checkCutoffForRestart, &feasible); 

[1315]  2732  if ((specialOptions_&262144) != 0 && !parentModel_) { 

 2733  // Save stuff and return! 

 2734  storedRowCuts_>saveStuff(bestObjective_, bestSolution_, 

 2735  solver_>getColLower(), 

 2736  solver_>getColUpper()); 

 2737  delete [] lowerBefore; 

 2738  delete [] upperBefore; 

 2739  delete saveSolver; 

 2740  return; 

 2741  } 

[1286]  2742  /* 

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

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

 2745  currently uses strong branching to evaluate branch object candidates, 

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

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

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

[2]  2749  

[1286]  2750  Monotone variables trigger a call to resolve(). If the problem remains 

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

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

[2]  2753  

[1286]  2754  Loss of feasibility will result in the deletion of newNode. 

 2755  */ 

[2]  2756  

[1286]  2757  bool resolved = false ; 

 2758  CbcNode *newNode = NULL ; 

 2759  numberFixedAtRoot_ = 0; 

 2760  numberFixedNow_ = 0; 

 2761  int numberIterationsAtContinuous = numberIterations_; 

 2762  //solverCharacteristics_>setSolver(solver_); 

 2763  if (feasible) { 

 2764  //#define HOTSTART 1 

[1271]  2765  #if HOTSTART<0 

[1286]  2766  if (bestSolution_ && !parentModel_ && !hotstartSolution_ && 

 2767  (moreSpecialOptions_&1024) != 0) { 

 2768  // Set priorities so only branch on ones we need to 

 2769  // use djs and see if only few branches needed 

[1271]  2770  #ifndef NDEBUG 

[1286]  2771  double integerTolerance = getIntegerTolerance() ; 

[1271]  2772  #endif 

[1286]  2773  bool possible = true; 

 2774  const double * saveLower = continuousSolver_>getColLower(); 

 2775  const double * saveUpper = continuousSolver_>getColUpper(); 

 2776  for (int i = 0; i < numberObjects_; i++) { 

 2777  const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]); 

 2778  if (thisOne) { 

 2779  int iColumn = thisOne>columnNumber(); 

 2780  if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) { 

 2781  possible = false; 

 2782  break; 

 2783  } 

 2784  } else { 

 2785  possible = false; 

 2786  break; 

 2787  } 

 2788  } 

 2789  if (possible) { 

 2790  OsiSolverInterface * solver = continuousSolver_>clone(); 

 2791  int numberColumns = solver>getNumCols(); 

 2792  for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) { 

 2793  double value = bestSolution_[iColumn] ; 

 2794  value = CoinMax(value, saveLower[iColumn]) ; 

 2795  value = CoinMin(value, saveUpper[iColumn]) ; 

 2796  value = floor(value + 0.5); 

 2797  if (solver>isInteger(iColumn)) { 

 2798  solver>setColLower(iColumn, value); 

 2799  solver>setColUpper(iColumn, value); 

 2800  } 

 2801  } 

 2802  solver>setHintParam(OsiDoDualInResolve, false, OsiHintTry); 

 2803  // objlim and all slack 

 2804  double direction = solver>getObjSense(); 

 2805  solver>setDblParam(OsiDualObjectiveLimit, 1.0e50*direction); 

 2806  CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver>getEmptyWarmStart()); 

 2807  solver>setWarmStart(basis); 

 2808  delete basis; 

 2809  bool changed = true; 

 2810  hotstartPriorities_ = new int [numberColumns]; 

 2811  for (int iColumn = 0; iColumn < numberColumns; iColumn++) 

 2812  hotstartPriorities_[iColumn] = 1; 

 2813  while (changed) { 

 2814  changed = false; 

 2815  solver>resolve(); 

 2816  if (!solver>isProvenOptimal()) { 

 2817  possible = false; 

 2818  break; 

 2819  } 

 2820  const double * dj = solver>getReducedCost(); 

 2821  const double * colLower = solver>getColLower(); 

 2822  const double * colUpper = solver>getColUpper(); 

 2823  const double * solution = solver>getColSolution(); 

 2824  int nAtLbNatural = 0; 

 2825  int nAtUbNatural = 0; 

 2826  int nZeroDj = 0; 

 2827  int nForced = 0; 

 2828  for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) { 

 2829  double value = solution[iColumn] ; 

 2830  value = CoinMax(value, saveLower[iColumn]) ; 

 2831  value = CoinMin(value, saveUpper[iColumn]) ; 

 2832  if (solver>isInteger(iColumn)) { 

 2833  assert(fabs(value  solution[iColumn]) <= integerTolerance) ; 

 2834  if (hotstartPriorities_[iColumn] == 1) { 

 2835  if (dj[iColumn] < 1.0e6) { 

 2836  // negative dj 

 2837  if (saveUpper[iColumn] == colUpper[iColumn]) { 

 2838  nAtUbNatural++; 

 2839  hotstartPriorities_[iColumn] = 2; 

 2840  solver>setColLower(iColumn, saveLower[iColumn]); 

 2841  solver>setColUpper(iColumn, saveUpper[iColumn]); 

 2842  } else { 

 2843  nForced++; 

 2844  } 

 2845  } else if (dj[iColumn] > 1.0e6) { 

 2846  // positive dj 

 2847  if (saveLower[iColumn] == colLower[iColumn]) { 

 2848  nAtLbNatural++; 

 2849  hotstartPriorities_[iColumn] = 2; 

 2850  solver>setColLower(iColumn, saveLower[iColumn]); 

 2851  solver>setColUpper(iColumn, saveUpper[iColumn]); 

 2852  } else { 

 2853  nForced++; 

 2854  } 

 2855  } else { 

 2856  // zero dj 

 2857  nZeroDj++; 

 2858  } 

 2859  } 

 2860  } 

 2861  } 

[1271]  2862  #ifdef CLP_INVESTIGATE 

[1286]  2863  printf("%d forced, %d naturally at lower, %d at upper  %d zero dj\n", 

 2864  nForced, nAtLbNatural, nAtUbNatural, nZeroDj); 

[1271]  2865  #endif 

[1286]  2866  if (nAtLbNatural  nAtUbNatural) { 

 2867  changed = true; 

 2868  } else { 

 2869  if (nForced + nZeroDj > 50  

 2870  (nForced + nZeroDj)*4 > numberIntegers_) 

 2871  possible = false; 

 2872  } 

 2873  } 

 2874  delete solver; 

 2875  } 

 2876  if (possible) { 

 2877  setHotstartSolution(bestSolution_); 

 2878  if (!saveCompare) { 

 2879  // create depth first comparison 

 2880  saveCompare = nodeCompare_; 

 2881  // depth first 

 2882  nodeCompare_ = new CbcCompareDepth(); 

 2883  tree_>setComparison(*nodeCompare_) ; 

 2884  } 

 2885  } else { 

 2886  delete [] hotstartPriorities_; 

 2887  hotstartPriorities_ = NULL; 

 2888  } 

 2889  } 

[1271]  2890  #endif 

 2891  #if HOTSTART>0 

[1286]  2892  if (hotstartSolution_ && !hotstartPriorities_) { 

 2893  // Set up hot start 

 2894  OsiSolverInterface * solver = solver_>clone(); 

 2895  double direction = solver_>getObjSense() ; 

 2896  int numberColumns = solver>getNumCols(); 

 2897  double * saveLower = CoinCopyOfArray(solver>getColLower(), numberColumns); 

 2898  double * saveUpper = CoinCopyOfArray(solver>getColUpper(), numberColumns); 

 2899  // move solution 

 2900  solver>setColSolution(hotstartSolution_); 

 2901  // point to useful information 

 2902  const double * saveSolution = testSolution_; 

 2903  testSolution_ = solver>getColSolution(); 

 2904  OsiBranchingInformation usefulInfo = usefulInformation(); 

 2905  testSolution_ = saveSolution; 

 2906  /* 

 2907  Run through the objects and use feasibleRegion() to set variable bounds 

 2908  so as to fix the variables specified in the objects at their value in this 

 2909  solution. Since the object list contains (at least) one object for every 

 2910  integer variable, this has the effect of fixing all integer variables. 

 2911  */ 

 2912  for (int i = 0; i < numberObjects_; i++) 

 2913  object_[i]>feasibleRegion(solver, &usefulInfo); 

 2914  solver>resolve(); 

 2915  assert (solver>isProvenOptimal()); 

 2916  double gap = CoinMax((solver>getObjValue()  solver_>getObjValue()) * direction, 0.0) ; 

 2917  const double * dj = solver>getReducedCost(); 

 2918  const double * colLower = solver>getColLower(); 

 2919  const double * colUpper = solver>getColUpper(); 

 2920  const double * solution = solver>getColSolution(); 

 2921  int nAtLbNatural = 0; 

 2922  int nAtUbNatural = 0; 

 2923  int nAtLbNaturalZero = 0; 

 2924  int nAtUbNaturalZero = 0; 

 2925  int nAtLbFixed = 0; 

 2926  int nAtUbFixed = 0; 

 2927  int nAtOther = 0; 

 2928  int nAtOtherNatural = 0; 

 2929  int nNotNeeded = 0; 

 2930  delete [] hotstartSolution_; 

 2931  hotstartSolution_ = new double [numberColumns]; 

 2932  delete [] hotstartPriorities_; 

 2933  hotstartPriorities_ = new int [numberColumns]; 

 2934  int * order = (int *) saveUpper; 

 2935  int nFix = 0; 

 2936  double bestRatio = COIN_DBL_MAX; 

 2937  for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) { 

 2938  double value = solution[iColumn] ; 

 2939  value = CoinMax(value, saveLower[iColumn]) ; 

 2940  value = CoinMin(value, saveUpper[iColumn]) ; 

 2941  double sortValue = COIN_DBL_MAX; 

 2942  if (solver>isInteger(iColumn)) { 

 2943  assert(fabs(value  solution[iColumn]) <= 1.0e5) ; 

 2944  double value2 = floor(value + 0.5); 

 2945  if (dj[iColumn] < 1.0e6) { 

 2946  // negative dj 

 2947  //assert (value2==colUpper[iColumn]); 

 2948  if (saveUpper[iColumn] == colUpper[iColumn]) { 

 2949  nAtUbNatural++; 

 2950  sortValue = 0.0; 

 2951  double value = dj[iColumn]; 

 2952  if (value > gap) 

 2953  nFix++; 

 2954  else if (gap < value*bestRatio) 

 2955  bestRatio = gap / value; 

 2956  if (saveLower[iColumn] != colLower[iColumn]) { 

 2957  nNotNeeded++; 

 2958  sortValue = 1.0e20; 

 2959  } 

 2960  } else if (saveLower[iColumn] == colUpper[iColumn]) { 

 2961  nAtLbFixed++; 

 2962  sortValue = dj[iColumn]; 

 2963  } else { 

 2964  nAtOther++; 

 2965  sortValue = 0.0; 

 2966  if (saveLower[iColumn] != colLower[iColumn] && 

 2967  saveUpper[iColumn] != colUpper[iColumn]) { 

 2968  nNotNeeded++; 

 2969  sortValue = 1.0e20; 

 2970  } 

 2971  } 

 2972  } else if (dj[iColumn] > 1.0e6) { 

 2973  // positive dj 

 2974  //assert (value2==colLower[iColumn]); 

 2975  if (saveLower[iColumn] == colLower[iColumn]) { 

 2976  nAtLbNatural++; 

 2977  sortValue = 0.0; 

 2978  double value = dj[iColumn]; 

 2979  if (value > gap) 

 2980  nFix++; 

 2981  else if (gap < value*bestRatio) 

 2982  bestRatio = gap / value; 

 2983  if (saveUpper[iColumn] != colUpper[iColumn]) { 

 2984  nNotNeeded++; 

 2985  sortValue = 1.0e20; 

 2986  } 

 2987  } else if (saveUpper[iColumn] == colLower[iColumn]) { 

 2988  nAtUbFixed++; 

 2989  sortValue = dj[iColumn]; 

 2990  } else { 

 2991  nAtOther++; 

 2992  sortValue = 0.0; 

 2993  if (saveLower[iColumn] != colLower[iColumn] && 

 2994  saveUpper[iColumn] != colUpper[iColumn]) { 

 2995  nNotNeeded++; 

 2996  sortValue = 1.0e20; 

 2997  } 

 2998  } 

 2999  } else { 

 3000  // zero dj 

 3001  if (value2 == saveUpper[iColumn]) { 

 3002  nAtUbNaturalZero++; 

 3003  sortValue = 0.0; 

 3004  if (saveLower[iColumn] != colLower[iColumn]) { 

 3005  nNotNeeded++; 

 3006  sortValue = 1.0e20; 

 3007  } 

 3008  } else if (value2 == saveLower[iColumn]) { 

 3009  nAtLbNaturalZero++; 

 3010  sortValue = 0.0; 

 3011  } else { 

 3012  nAtOtherNatural++; 

 3013  sortValue = 0.0; 

 3014  if (saveLower[iColumn] != colLower[iColumn] && 

 3015  saveUpper[iColumn] != colUpper[iColumn]) { 

 3016  nNotNeeded++; 

 3017  sortValue = 1.0e20; 

 3018  } 

 3019  } 

 3020  } 

[1016]  3021  #if HOTSTART==3 

[1286]  3022  sortValue = fabs(dj[iColumn]); 

[1016]  3023  #endif 

[1286]  3024  } 

 3025  hotstartSolution_[iColumn] = value ; 

 3026  saveLower[iColumn] = sortValue; 

 3027  order[iColumn] = iColumn; 

 3028  } 

 3029  printf("** can fix %d columns  best ratio for others is %g on gap of %g\n", 

 3030  nFix, bestRatio, gap); 

 3031  int nNeg = 0; 

 3032  CoinSort_2(saveLower, saveLower + numberColumns, order); 

 3033  for (int i = 0; i < numberColumns; i++) { 

 3034  if (saveLower[i] < 0.0) { 

 3035  nNeg++; 

[1016]  3036  #if HOTSTART==2HOTSTART==3 

[1286]  3037  // swap sign ? 

 3038  saveLower[i] = saveLower[i]; 

[1016]  3039  #endif 

[1286]  3040  } 

 3041  } 

 3042  CoinSort_2(saveLower, saveLower + nNeg, order); 

 3043  for (int i = 0; i < numberColumns; i++) { 

[1016]  3044  #if HOTSTART==1 

[1286]  3045  hotstartPriorities_[order[i]] = 100; 

[1016]  3046  #else 

[1286]  3047  hotstartPriorities_[order[i]] = (i + 1); 

[1016]  3048  #endif 

[1286]  3049  } 

 3050  printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n", 

 3051  nAtLbNatural, 

 3052  nAtUbNatural, 

 3053  nAtLbNaturalZero, 

 3054  nAtUbNaturalZero, 

 3055  nAtLbFixed, 

 3056  nAtUbFixed, 

 3057  nAtOther, 

 3058  nAtOtherNatural, nNotNeeded, nNeg); 

 3059  delete [] saveLower; 

 3060  delete [] saveUpper; 

 3061  if (!saveCompare) { 

 3062  // create depth first comparison 

 3063  saveCompare = nodeCompare_; 

 3064  // depth first 

 3065  nodeCompare_ = new CbcCompareDepth(); 

 3066  tree_>setComparison(*nodeCompare_) ; 

 3067  } 

 3068  } 

[1016]  3069  #endif 

[1286]  3070  newNode = new CbcNode ; 

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

 3072  setObjectiveValue(newNode, NULL); 

 3073  anyAction = 1 ; 

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

 3075  CbcNode fakeNode; 

 3076  if (!currentNode_) { 

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

 3078  currentNode_ = &fakeNode; 

 3079  } 

 3080  phase_ = 3; 

 3081  // only allow 1000 passes 

 3082  int numberPassesLeft = 1000; 

 3083  // This is first crude step 

 3084  if (numberAnalyzeIterations_) { 

 3085  delete [] analyzeResults_; 

 3086  analyzeResults_ = new double [4*numberIntegers_]; 

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

 3088  if (numberFixedAtRoot_ > 0) { 

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

 3090  setPointers(solver_); 

 3091  numberFixedNow_ = numberFixedAtRoot_; 

 3092  } else if (numberFixedAtRoot_ < 0) { 

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

 3094  anyAction = 2; 

 3095  delete newNode ; 

 3096  newNode = NULL ; 

 3097  feasible = false ; 

 3098  } 

 3099  } 

 3100  OsiSolverBranch * branches = NULL; 

 3101  anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved, 

 3102  NULL, NULL, NULL, branches); 

 3103  if (anyAction == 2  newNode>objectiveValue() >= cutoff) { 

 3104  if (anyAction != 2) { 

 3105  // zap parent nodeInfo 

[640]  3106  #ifdef COIN_DEVELOP 

[1286]  3107  printf("zapping CbcNodeInfo %x\n", reinterpret_cast<int>(newNode>nodeInfo()>parent())); 

[640]  3108  #endif 

[1286]  3109  if (newNode>nodeInfo()) 

 3110  newNode>nodeInfo()>nullParent(); 

 3111  } 

 3112  delete newNode ; 

 3113  newNode = NULL ; 

 3114  feasible = false ; 

 3115  } 

[640]  3116  } 

[1315]  3117  if (newNode && probingInfo_) { 

 3118  int number01 = probingInfo_>numberIntegers(); 

 3119  //const fixEntry * entry = probingInfo_>fixEntries(); 

 3120  const int * toZero = probingInfo_>toZero(); 

 3121  //const int * toOne = probingInfo_>toOne(); 

 3122  //const int * integerVariable = probingInfo_>integerVariable(); 

 3123  if (toZero[number01]) { 

 3124  CglTreeProbingInfo info(*probingInfo_); 

[1393]  3125  #ifdef JJF_ZERO 

[1409]  3126  /* 

 3127  Marginal idea. Further exploration probably good. Build some extra 

 3128  cliques from probing info. Not quite worth the effort? 

 3129  */ 

[1315]  3130  OsiSolverInterface * fake = info.analyze(*solver_, 1); 

 3131  if (fake) { 

 3132  fake>writeMps("fake"); 

 3133  CglFakeClique cliqueGen(fake); 

 3134  cliqueGen.setStarCliqueReport(false); 

 3135  cliqueGen.setRowCliqueReport(false); 

 3136  cliqueGen.setMinViolation(0.1); 

 3137  addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, 200); 

 3138  generator_[numberCutGenerators_1]>setTiming(true); 

 3139  } 

 3140  #endif 

 3141  if (probingInfo_>packDown()) { 

 3142  #ifdef CLP_INVESTIGATE 

 3143  printf("%d implications on %d 01\n", toZero[number01], number01); 

 3144  #endif 

[1409]  3145  // Create a cut generator that remembers implications discovered at root. 

[1315]  3146  CglImplication implication(probingInfo_); 

 3147  addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, 200); 

 3148  } else { 

 3149  delete probingInfo_; 

 3150  probingInfo_ = NULL; 

 3151  } 

 3152  } else { 

 3153  delete probingInfo_; 

 3154  

 3155  probingInfo_ = NULL; 

 3156  } 

 3157  } 

[1286]  3158  /* 

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

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

 3161  current objective cutoff. 

 3162  */ 

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

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

 3165  // initialize for print out 

 3166  int lastDepth = 0; 

 3167  int lastUnsatisfied = 0; 

 3168  if (newNode) 

 3169  lastUnsatisfied = newNode>numberUnsatisfied(); 

 3170  /* 

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

 3172  integrality (i.e., newNode>branchingObject(), indicating we've been able to 

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

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

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

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

 3177  (newNode>branchingObject() is NULL), takeOffCuts() will ensure that the solver holds 

 3178  a valid solution for use by setBestSolution(). 

 3179  */ 

 3180  CoinWarmStartBasis *lastws = NULL ; 

 3181  if (feasible && newNode>branchingObject()) { 

 3182  if (resolved) { 

 3183  takeOffCuts(cuts, false, NULL) ; 

[2]  3184  # ifdef CHECK_CUT_COUNTS 

[1286]  3185  { printf("Number of rows after chooseBranch fix (root)" 

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

 3187  numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ; 

 3188  const CoinWarmStartBasis* debugws = 

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

 3190  debugws>print() ; 

 3191  delete debugws ; 

 3192  } 

[2]  3193  # endif 

[1286]  3194  } 

 3195  //newNode>createInfo(this,NULL,NULL,NULL,NULL,0,0) ; 

 3196  //newNode>nodeInfo()>addCuts(cuts, 

 3197  // newNode>numberBranches(),whichGenerator_) ; 

 3198  if (lastws) delete lastws ; 

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

[2]  3200  } 

[1286]  3201  /* 

 3202  Continuous data to be used later 

 3203  */ 

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

 3205  continuousInfeasibilities_ = 0 ; 

 3206  if (newNode) { 

 3207  continuousObjective_ = newNode>objectiveValue() ; 

 3208  delete [] continuousSolution_; 

 3209  continuousSolution_ = CoinCopyOfArray(solver_>getColSolution(), 

 3210  numberColumns); 

 3211  continuousInfeasibilities_ = newNode>numberUnsatisfied() ; 

 3212  } 

 3213  /* 

 3214  Bound may have changed so reset in objects 

 3215  */ 

 3216  { 

 3217  int i ; 

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

 3219  object_[i]>resetBounds(solver_) ; 

 3220  } 

 3221  /* 

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

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

 3224  integral solution at the root node. 

[2]  3225  

[1286]  3226  initializeInfo sets the reference counts in the nodeInfo object. Since 

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

 3228  */ 

 3229  double bestValue = 0.0 ; 

 3230  if (newNode) { 

 3231  bestValue = newNode>objectiveValue(); 

 3232  if (newNode>branchingObject()) { 

 3233  newNode>initializeInfo() ; 

 3234  tree_>push(newNode) ; 

 3235  if (statistics_) { 

 3236  if (numberNodes2_ == maximumStatistics_) { 

 3237  maximumStatistics_ = 2 * maximumStatistics_; 

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

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

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

 3241  delete [] statistics_; 

 3242  statistics_ = temp; 

 3243  } 

 3244  assert (!statistics_[numberNodes2_]); 

 3245  statistics_[numberNodes2_] = new CbcStatistics(newNode, this); 

 3246  } 

 3247  numberNodes2_++; 

[2]  3248  # ifdef CHECK_NODE 

[1286]  3249  printf("Node %x on tree\n", newNode) ; 

[2]  3250  # endif 

[1286]  3251  } else { 

 3252  // continuous is integer 

 3253  double objectiveValue = newNode>objectiveValue(); 

 3254  setBestSolution(CBC_SOLUTION, objectiveValue, 

 3255  solver_>getColSolution()) ; 

 3256  delete newNode ; 

 3257  newNode = NULL ; 

 3258  } 

[2]  3259  } 

 3260  

[1286]  3261  if (printFrequency_ <= 0) { 

 3262  printFrequency_ = 1000 ; 

 3263  if (getNumCols() > 2000) 

 3264  printFrequency_ = 100 ; 

 3265  } 

 3266  /* 

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

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

 3269  */ 

 3270  numberLongStrong_ = 0; 

 3271  CbcNode * createdNode = NULL; 

[1122]  3272  #ifdef CBC_THREAD 

[1409]  3273  if ((specialOptions_&2048) != 0) 

 3274  numberThreads_ = 0; 

 3275  if (numberThreads_ ) { 

[1286]  3276  nodeCompare_>sayThreaded(); // need to use addresses 

[1409]  3277  master_ = new CbcBaseModel(*this, 

 3278  (parallelMode() < 1) ? 1 : 0); 

 3279  masterThread_ = master_>masterThread(); 

[687]  3280  } 

 3281  #endif 

[931]  3282  #ifdef COIN_HAS_CLP 

[1286]  3283  { 

 3284  OsiClpSolverInterface * clpSolver 

 3285  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 3286  if (clpSolver && !parentModel_) { 

 3287  clpSolver>computeLargestAway(); 

 3288  } 

[931]  3289  } 

 3290  #endif 

[1286]  3291  /* 

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

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

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

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

 3296  add the node to the live set. 

[2]  3297  

[1286]  3298  The first action is to winnow the live set to remove nodes which are worse 

 3299  than the current objective cutoff. 

 3300  */ 

 3301  if (solver_>getRowCutDebuggerAlways()) { 

 3302  OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_>getRowCutDebuggerAlways()); 

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

 3304  if (!debugger) { 

 3305  // infeasible!! 

 3306  printf("before search\n"); 

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

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

 3309  const double * solution = debuggerX>optimalSolution(); 

 3310  int numberColumns = solver_>getNumCols(); 

 3311  for (int i = 0; i < numberColumns; i++) { 

 3312  if (solver_>isInteger(i)) { 

 3313  if (solution[i] < lower[i]  1.0e6  solution[i] > upper[i] + 1.0e6) 

 3314  printf("**** "); 

 3315  printf("%d %g <= %g <= %g\n", 

 3316  i, lower[i], solution[i], upper[i]); 

 3317  } 

 3318  } 

 3319  //abort(); 

 3320  } 

[789]  3321  } 

[1286]  3322  { 

 3323  // may be able to change cutoff now 

 3324  double cutoff = getCutoff(); 

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

 3326  if (cutoff > bestObjective_  increment) { 

 3327  cutoff = bestObjective_  increment ; 

 3328  setCutoff(cutoff) ; 

 3329  } 

[1271]  3330  } 

[1121]  3331  #ifdef CBC_THREAD 

[1286]  3332  bool goneParallel = false; 

[1121]  3333  #endif 

[838]  3334  #define MAX_DEL_NODE 1 

[1286]  3335  CbcNode * delNode[MAX_DEL_NODE+1]; 

 3336  int nDeleteNode = 0; 

 3337  // For Printing etc when parallel 

 3338  int lastEvery1000 = 0; 

 3339  int lastPrintEvery = 0; 

[1315]  3340  int numberConsecutiveInfeasible = 0; 

[1286]  3341  while (true) { 

[1409]  3342  lockThread(); 

[931]  3343  #ifdef COIN_HAS_CLP 

[1409]  3344  // See if we want dantzig row choice 

 3345  goToDantzig(100, savePivotMethod); 

[1357]  3346  #endif 

[1286]  3347  if (tree_>empty()) { 

[1122]  3348  #ifdef CBC_THREAD 

[1412]  3349  if (parallelMode() > 0 && master_) { 

[1409]  3350  int anyLeft = master_>waitForThreadsInTree(0); 

[1412]  3351  if (!anyLeft) { 

 3352  master_>stopThreads(1); 

[1409]  3353  break; 

[1412]  3354  } 

[1409]  3355  } else { 

 3356  break; 

[1286]  3357  } 

[1409]  3358  #else 

 3359  break; 

[1122]  3360  #endif 

[1409]  3361  } else { 

[1286]  3362  unlockThread(); 

 3363  } 

 3364  // If done 100 nodes see if worth trying reduction 

 3365  if (numberNodes_ == 50  numberNodes_ == 100) { 

[1271]  3366  #ifdef COIN_HAS_CLP 

[1286]  3367  OsiClpSolverInterface * clpSolver 

 3368  = dynamic_cast<OsiClpSolverInterface *> (solver_); 

 3369  if (clpSolver && ((specialOptions_&131072) == 0) && true) { 

 3370  ClpSimplex * simplex = clpSolver>getModelPtr(); 

 3371  int perturbation = simplex>perturbation(); 

[1132]  3372  #ifdef CLP_INVESTIGATE 

[1286]  3373  printf("Testing its n,s %d %d solves n,s %d %d  pert %d\n", 

 3374  numberIterations_, saveNumberIterations, 

 3375  numberSolves_, saveNumberSolves, perturbation); 

[1132]  3376  #endif 

[1286]  3377  if (perturbation == 50 && (numberIterations_  saveNumberIterations) < 

 3378  8*(numberSolves_  saveNumberSolves)) { 

 3379  // switch off perturbation 

 3380  simplex>setPerturbation(100); 

[1271]  3381  #ifdef PERTURB_IN_FATHOM 

[1286]  3382  // but allow in fathom 

 3383  specialOptions_ = 131072; 

[1271]  3384  #endif 

 3385  #ifdef CLP_INVESTIGATE 

[1286]  3386  printf("Perturbation switched off\n"); 

[1271]  3387  #endif 

[1286]  3388  } 

 3389  } 


