source: trunk/Cbc/src/CbcModel.cpp @ 2464

Last change on this file since 2464 was 2464, checked in by unxusr, 2 years ago

.clang-format with proposal for formatting code

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 702.4 KB
Line 
1/* $Id: CbcModel.cpp 2464 2019-01-03 19:03:23Z unxusr $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#pragma warning(disable : 4786)
9#endif
10
11#include "CbcConfig.h"
12
13#include <string>
14//#define CBC_DEBUG 1
15//#define CHECK_CUT_COUNTS
16//#define CHECK_NODE
17//#define CHECK_NODE_FULL
18//#define NODE_LOG
19//#define GLOBAL_CUTS_JUST_POINTERS
20#ifdef CGL_DEBUG_GOMORY
21extern int gomory_try;
22#endif
23#include <cassert>
24#include <cmath>
25#include <cfloat>
26#ifdef COIN_HAS_CLP
27// include Presolve from Clp
28#include "ClpPresolve.hpp"
29#include "OsiClpSolverInterface.hpp"
30#include "ClpNode.hpp"
31#include "ClpDualRowDantzig.hpp"
32#include "ClpSimplexPrimal.hpp"
33#endif
34
35#include "CbcEventHandler.hpp"
36
37#include "OsiSolverInterface.hpp"
38#include "OsiAuxInfo.hpp"
39#include "OsiSolverBranch.hpp"
40#include "OsiChooseVariable.hpp"
41#include "CoinWarmStartBasis.hpp"
42#include "CoinPackedMatrix.hpp"
43#include "CoinHelperFunctions.hpp"
44#include "CbcBranchActual.hpp"
45#include "CbcBranchDynamic.hpp"
46#include "CbcHeuristic.hpp"
47#include "CbcHeuristicFPump.hpp"
48#include "CbcHeuristicRINS.hpp"
49#include "CbcHeuristicDive.hpp"
50#include "CbcModel.hpp"
51#include "CbcTreeLocal.hpp"
52#include "CbcStatistics.hpp"
53#include "CbcStrategy.hpp"
54#include "CbcMessage.hpp"
55#include "OsiRowCut.hpp"
56#include "OsiColCut.hpp"
57#include "OsiRowCutDebugger.hpp"
58#include "OsiCuts.hpp"
59#include "CbcCountRowCut.hpp"
60#include "CbcCutGenerator.hpp"
61#include "CbcFeasibilityBase.hpp"
62#include "CbcFathom.hpp"
63#include "CbcFullNodeInfo.hpp"
64#ifdef COIN_HAS_NTY
65#include "CbcSymmetry.hpp"
66#endif
67// include Probing
68#include "CglProbing.hpp"
69#include "CglGomory.hpp"
70#include "CglTwomir.hpp"
71// include preprocessing
72#include "CglPreProcess.hpp"
73#include "CglDuplicateRow.hpp"
74#include "CglStored.hpp"
75#include "CglClique.hpp"
76#include "CglKnapsackCover.hpp"
77
78#include "CoinTime.hpp"
79#include "CoinMpsIO.hpp"
80
81#include "CbcCompareActual.hpp"
82#include "CbcTree.hpp"
83// This may be dummy
84#include "CbcThread.hpp"
85/* Various functions local to CbcModel.cpp */
86
87typedef struct {
88  double useCutoff;
89  CbcModel *model;
90  int switches;
91} rootBundle;
92static void *doRootCbcThread(void *voidInfo);
93
94namespace {
95
96//-------------------------------------------------------------------
97// Returns the greatest common denominator of two
98// positive integers, a and b, found using Euclid's algorithm
99//-------------------------------------------------------------------
100static int gcd(int a, int b)
101{
102  int remainder = -1;
103  // make sure a<=b (will always remain so)
104  if (a > b) {
105    // Swap a and b
106    int temp = a;
107    a = b;
108    b = temp;
109  }
110  // if zero then gcd is nonzero (zero may occur in rhs of packed)
111  if (!a) {
112    if (b) {
113      return b;
114    } else {
115      printf("**** gcd given two zeros!!\n");
116      abort();
117    }
118  }
119  while (remainder) {
120    remainder = b % a;
121    b = a;
122    a = remainder;
123  }
124  return b;
125}
126
127#ifdef CHECK_NODE_FULL
128
129/*
130  Routine to verify that tree linkage is correct. The invariant that is tested
131  is
132
133  reference count = (number of actual references) + (number of branches left)
134
135  The routine builds a set of paired arrays, info and count, by traversing the
136  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
137  referenced (via the parent field) is recorded in count. Then a final check is
138  made to see if the numberPointingToThis_ field agrees.
139*/
140
141void verifyTreeNodes(const CbcTree *branchingTree, const CbcModel &model)
142
143{
144  if (model.getNodeCount() == 661)
145    return;
146  printf("*** CHECKING tree after %d nodes\n", model.getNodeCount());
147
148  int j;
149  int nNodes = branchingTree->size();
150#define MAXINFO 1000
151  int *count = new int[MAXINFO];
152  CbcNodeInfo **info = new CbcNodeInfo *[MAXINFO];
153  int nInfo = 0;
154  /*
155      Collect all CbcNodeInfo objects in info, by starting from each live node and
156      traversing back to the root. Nodes in the live set should have unexplored
157      branches remaining.
158
159      TODO: The `while (nodeInfo)' loop could be made to break on reaching a
160        common ancester (nodeInfo is found in info[k]). Alternatively, the
161        check could change to signal an error if nodeInfo is not found above a
162        common ancestor.
163    */
164  for (j = 0; j < nNodes; j++) {
165    CbcNode *node = branchingTree->nodePointer(j);
166    if (!node)
167      continue;
168    CbcNodeInfo *nodeInfo = node->nodeInfo();
169    int change = node->nodeInfo()->numberBranchesLeft();
170    assert(change);
171    while (nodeInfo) {
172      int k;
173      for (k = 0; k < nInfo; k++) {
174        if (nodeInfo == info[k])
175          break;
176      }
177      if (k == nInfo) {
178        assert(nInfo < MAXINFO);
179        nInfo++;
180        info[k] = nodeInfo;
181        count[k] = 0;
182      }
183      nodeInfo = nodeInfo->parent();
184    }
185  }
186  /*
187      Walk the info array. For each nodeInfo, look up its parent in info and
188      increment the corresponding count.
189    */
190  for (j = 0; j < nInfo; j++) {
191    CbcNodeInfo *nodeInfo = info[j];
192    nodeInfo = nodeInfo->parent();
193    if (nodeInfo) {
194      int k;
195      for (k = 0; k < nInfo; k++) {
196        if (nodeInfo == info[k])
197          break;
198      }
199      assert(k < nInfo);
200      count[k]++;
201    }
202  }
203  /*
204      Walk the info array one more time and check that the invariant holds. The
205      number of references (numberPointingToThis()) should equal the sum of the
206      number of actual references (held in count[]) plus the number of potential
207      references (unexplored branches, numberBranchesLeft()).
208    */
209  for (j = 0; j < nInfo; j++) {
210    CbcNodeInfo *nodeInfo = info[j];
211    if (nodeInfo) {
212      int k;
213      for (k = 0; k < nInfo; k++)
214        if (nodeInfo == info[k])
215          break;
216      printf("Nodeinfo %x - %d left, %d count\n",
217        nodeInfo,
218        nodeInfo->numberBranchesLeft(),
219        nodeInfo->numberPointingToThis());
220      assert(nodeInfo->numberPointingToThis() == count[k] + nodeInfo->numberBranchesLeft());
221    }
222  }
223
224  delete[] count;
225  delete[] info;
226
227  return;
228}
229
230#endif /* CHECK_NODE_FULL */
231
232#ifdef CHECK_CUT_COUNTS
233
234/*
235  Routine to verify that cut reference counts are correct.
236*/
237void verifyCutCounts(const CbcTree *branchingTree, CbcModel &model)
238
239{
240  printf("*** CHECKING cuts after %d nodes\n", model.getNodeCount());
241
242  int j;
243  int nNodes = branchingTree->size();
244
245  /*
246      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
247      in all cuts. We traverse the tree by starting from each live node and working
248      back to the root. At each CbcNodeInfo, check for cuts.
249    */
250  for (j = 0; j < nNodes; j++) {
251    CbcNode *node = branchingTree->nodePointer(j);
252    CbcNodeInfo *nodeInfo = node->nodeInfo();
253    assert(node->nodeInfo()->numberBranchesLeft());
254    while (nodeInfo) {
255      int k;
256      for (k = 0; k < nodeInfo->numberCuts(); k++) {
257        CbcCountRowCut *cut = nodeInfo->cuts()[k];
258        if (cut)
259          cut->tempNumber_ = 0;
260      }
261      nodeInfo = nodeInfo->parent();
262    }
263  }
264  /*
265      Walk the live set again, this time collecting the list of cuts in use at each
266      node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
267      that when we recreate the basis for a node, we compress out the slack cuts.
268    */
269  for (j = 0; j < nNodes; j++) {
270    CoinWarmStartBasis *debugws = model.getEmptyBasis();
271    CbcNode *node = branchingTree->nodePointer(j);
272    CbcNodeInfo *nodeInfo = node->nodeInfo();
273    int change = node->nodeInfo()->numberBranchesLeft();
274    printf("Node %d %x (info %x) var %d way %d obj %g", j, node,
275      node->nodeInfo(), node->columnNumber(), node->way(),
276      node->objectiveValue());
277
278    model.addCuts1(node, debugws);
279
280    int i;
281    int numberRowsAtContinuous = model.numberRowsAtContinuous();
282    CbcCountRowCut **addedCuts = model.addedCuts();
283    for (i = 0; i < model.currentNumberCuts(); i++) {
284      CoinWarmStartBasis::Status status = debugws->getArtifStatus(i + numberRowsAtContinuous);
285      if (status != CoinWarmStartBasis::basic && addedCuts[i]) {
286        addedCuts[i]->tempNumber_ += change;
287      }
288    }
289
290    while (nodeInfo) {
291      nodeInfo = nodeInfo->parent();
292      if (nodeInfo)
293        printf(" -> %x", nodeInfo);
294    }
295    printf("\n");
296    delete debugws;
297  }
298  /*
299      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.
300
301      TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
302    */
303  for (j = 0; j < nNodes; j++) {
304    CbcNode *node = branchingTree->nodePointer(j);
305    CbcNodeInfo *nodeInfo = node->nodeInfo();
306    while (nodeInfo) {
307      int k;
308      for (k = 0; k < nodeInfo->numberCuts(); k++) {
309        CbcCountRowCut *cut = nodeInfo->cuts()[k];
310        if (cut && cut->tempNumber_ >= 0) {
311          if (cut->tempNumber_ != cut->numberPointingToThis())
312            printf("mismatch %x %d %x %d %d\n", nodeInfo, k,
313              cut, cut->tempNumber_, cut->numberPointingToThis());
314          else
315            printf("   match %x %d %x %d %d\n", nodeInfo, k,
316              cut, cut->tempNumber_, cut->numberPointingToThis());
317          cut->tempNumber_ = -1;
318        }
319      }
320      nodeInfo = nodeInfo->parent();
321    }
322  }
323
324  return;
325}
326
327#endif /* CHECK_CUT_COUNTS */
328
329#ifdef CHECK_CUT_SIZE
330
331/*
332  Routine to verify that cut reference counts are correct.
333*/
334void verifyCutSize(const CbcTree *branchingTree, CbcModel &model)
335{
336
337  int j;
338  int nNodes = branchingTree->size();
339  int totalCuts = 0;
340
341  /*
342      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
343      in all cuts. We traverse the tree by starting from each live node and working
344      back to the root. At each CbcNodeInfo, check for cuts.
345    */
346  for (j = 0; j < nNodes; j++) {
347    CbcNode *node = branchingTree->nodePointer(j);
348    CbcNodeInfo *nodeInfo = node->nodeInfo();
349    assert(node->nodeInfo()->numberBranchesLeft());
350    while (nodeInfo) {
351      totalCuts += nodeInfo->numberCuts();
352      nodeInfo = nodeInfo->parent();
353    }
354  }
355  printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n", model.getNodeCount(), totalCuts);
356  return;
357}
358
359#endif /* CHECK_CUT_SIZE */
360
361}
362
363/* End unnamed namespace for CbcModel.cpp */
364
365void CbcModel::analyzeObjective()
366/*
367  Try to find a minimum change in the objective function. The first scan
368  checks that there are no continuous variables with non-zero coefficients,
369  and grabs the largest objective coefficient associated with an unfixed
370  integer variable. The second scan attempts to scale up the objective
371  coefficients to a point where they are sufficiently close to integer that
372  we can pretend they are integer, and calculate a gcd over the coefficients
373  of interest. This will be the minimum increment for the scaled coefficients.
374  The final action is to scale the increment back for the original coefficients
375  and install it, if it's better than the existing value.
376
377  John's note: We could do better than this.
378
379  John's second note - apologies for changing s to z
380*/
381{
382  const double *objective = getObjCoefficients();
383  const double *lower = getColLower();
384  const double *upper = getColUpper();
385  /*
386      Scan continuous and integer variables to see if continuous
387      are cover or network with integral rhs.
388    */
389  double continuousMultiplier = 1.0;
390  double *coeffMultiplier = NULL;
391  double largestObj = 0.0;
392  double smallestObj = COIN_DBL_MAX;
393  {
394    const double *rowLower = getRowLower();
395    const double *rowUpper = getRowUpper();
396    int numberRows = solver_->getNumRows();
397    double *rhs = new double[numberRows];
398    memset(rhs, 0, numberRows * sizeof(double));
399    int iColumn;
400    int numberColumns = solver_->getNumCols();
401    // Column copy of matrix
402    int problemType = -1;
403    const double *element = solver_->getMatrixByCol()->getElements();
404    const int *row = solver_->getMatrixByCol()->getIndices();
405    const CoinBigIndex *columnStart = solver_->getMatrixByCol()->getVectorStarts();
406    const int *columnLength = solver_->getMatrixByCol()->getVectorLengths();
407    int numberInteger = 0;
408    int numberIntegerObj = 0;
409    int numberGeneralIntegerObj = 0;
410    int numberIntegerWeight = 0;
411    int numberContinuousObj = 0;
412    double cost = COIN_DBL_MAX;
413    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
414      if (upper[iColumn] == lower[iColumn]) {
415        CoinBigIndex start = columnStart[iColumn];
416        CoinBigIndex end = start + columnLength[iColumn];
417        for (CoinBigIndex j = start; j < end; j++) {
418          int iRow = row[j];
419          rhs[iRow] += lower[iColumn] * element[j];
420        }
421      } else {
422        double objValue = objective[iColumn];
423        if (solver_->isInteger(iColumn))
424          numberInteger++;
425        if (objValue) {
426          if (!solver_->isInteger(iColumn)) {
427            numberContinuousObj++;
428          } else {
429            largestObj = CoinMax(largestObj, fabs(objValue));
430            smallestObj = CoinMin(smallestObj, fabs(objValue));
431            numberIntegerObj++;
432            if (cost == COIN_DBL_MAX)
433              cost = objValue;
434            else if (cost != objValue)
435              cost = -COIN_DBL_MAX;
436            int gap = static_cast<int>(upper[iColumn] - lower[iColumn]);
437            if (gap > 1) {
438              numberGeneralIntegerObj++;
439              numberIntegerWeight += gap;
440            }
441          }
442        }
443      }
444    }
445    int iType = 0;
446    if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 && numberIntegerObj * 3 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100)
447      iType = 1 + 4 + (((moreSpecialOptions_ & 536870912) == 0) ? 2 : 0);
448    else if (!numberContinuousObj && numberIntegerObj <= 100 && numberIntegerObj * 5 < numberObjects_ && numberIntegerWeight <= 100 && !parentModel_ && solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
449      iType = 4 + (((moreSpecialOptions_ & 536870912) == 0) ? 2 : 0);
450    else if (!numberContinuousObj && numberIntegerObj <= 100 && numberIntegerObj * 5 < numberObjects_ && !parentModel_ && 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<=5||cost!=-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        // no ? 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      }
517#ifdef COIN_HAS_CLP
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();
531#ifdef JJF_ZERO
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.0e-6 || value > upperValue - 1.0e-6)
541              nAtBound++;
542          }
543        }
544#endif
545        /*
546                  Generate a random objective function for problems where the given objective
547                  function is not terribly useful. (Nearly feasible, single integer variable,
548                  that sort of thing.
549                */
550        CoinDrand48(true, 1234567);
551        for (int i = 0; i < numberColumns; i++) {
552          double lowerValue = lower[i];
553          double upperValue = upper[i];
554          double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000)
555                                      : i + 1 + columnLength[i] * 1000;
556          value *= 0.001;
557          //value += columnLength[i];
558          if (lowerValue > -1.0e5 || upperValue < 1.0e5) {
559            if (fabs(lowerValue) > fabs(upperValue))
560              value = -value;
561            if (clpSolver->isInteger(i)) {
562              double solValue = solution[i];
563              // Better to add in 0.5 or 1.0??
564              if (solValue < lowerValue + 1.0e-6)
565                value = fabs(value) + 0.5; //fabs(value*1.5);
566              else if (solValue > upperValue - 1.0e-6)
567                value = -fabs(value) - 0.5; //-fabs(value*1.5);
568            }
569          } else {
570            value = 0.0;
571          }
572          fakeObj[i] = value;
573        }
574        // pass to solver
575        clpSolver->setFakeObjective(fakeObj);
576        delete[] fakeObj;
577      }
578#endif
579    } else if (largestObj < smallestObj * 5.0 && !parentModel_ && !numberContinuousObj && !numberGeneralIntegerObj && numberIntegerObj * 2 < CoinMin(numberColumns, 20)) {
580      // up priorities on costed
581      int iPriority = -1;
582      for (int i = 0; i < numberObjects_; i++) {
583        int k = object_[i]->priority();
584        if (iPriority == -1)
585          iPriority = k;
586        else if (iPriority != k)
587          iPriority = -2;
588      }
589      if (iPriority >= 100) {
590#if CBC_USEFUL_PRINTING > 1
591        printf("Setting variables with obj to high priority\n");
592#endif
593        for (int i = 0; i < numberObjects_; i++) {
594          CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(object_[i]);
595          if (obj) {
596            int iColumn = obj->columnNumber();
597            if (objective[iColumn])
598              object_[i]->setPriority(iPriority - 1);
599          }
600        }
601      }
602    }
603    int iRow;
604    for (iRow = 0; iRow < numberRows; iRow++) {
605      if (rowLower[iRow] > -1.0e20 && fabs(rowLower[iRow] - rhs[iRow] - floor(rowLower[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
606        continuousMultiplier = 0.0;
607        break;
608      }
609      if (rowUpper[iRow] < 1.0e20 && fabs(rowUpper[iRow] - rhs[iRow] - floor(rowUpper[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
610        continuousMultiplier = 0.0;
611        break;
612      }
613      // set rhs to limiting value
614      if (rowLower[iRow] != rowUpper[iRow]) {
615        if (rowLower[iRow] > -1.0e20) {
616          if (rowUpper[iRow] < 1.0e20) {
617            // no good
618            continuousMultiplier = 0.0;
619            break;
620          } else {
621            rhs[iRow] = rowLower[iRow] - rhs[iRow];
622            if (problemType < 0)
623              problemType = 3; // set cover
624            else if (problemType != 3)
625              problemType = 4;
626          }
627        } else {
628          rhs[iRow] = rowUpper[iRow] - rhs[iRow];
629          if (problemType < 0)
630            problemType = 1; // set partitioning <=
631          else if (problemType != 1)
632            problemType = 4;
633        }
634      } else {
635        rhs[iRow] = rowUpper[iRow] - rhs[iRow];
636        if (problemType < 0)
637          problemType = 3; // set partitioning ==
638        else if (problemType != 2)
639          problemType = 2;
640      }
641      if (fabs(rhs[iRow] - 1.0) > 1.0e-12)
642        problemType = 4;
643    }
644    if (continuousMultiplier) {
645      // 1 network, 2 cover, 4 negative cover
646      int possible = 7;
647      bool unitRhs = true;
648      // See which rows could be set cover
649      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
650        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
651          CoinBigIndex start = columnStart[iColumn];
652          CoinBigIndex end = start + columnLength[iColumn];
653          for (CoinBigIndex j = start; j < end; j++) {
654            double value = element[j];
655            if (value == 1.0) {
656            } else if (value == -1.0) {
657              rhs[row[j]] = -0.5;
658            } else {
659              rhs[row[j]] = -COIN_DBL_MAX;
660            }
661          }
662        }
663      }
664      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
665        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
666          if (!isInteger(iColumn)) {
667            CoinBigIndex start = columnStart[iColumn];
668            CoinBigIndex end = start + columnLength[iColumn];
669            double rhsValue = 0.0;
670            // 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
671            int type = 0;
672            for (CoinBigIndex j = start; j < end; j++) {
673              double value = element[j];
674              if (fabs(value) != 1.0) {
675                type = 3;
676                break;
677              } else if (value == 1.0) {
678                if (!type)
679                  type = 1;
680                else if (type != 1)
681                  type = 2;
682              } else {
683                if (!type)
684                  type = -1;
685                else if (type != -1)
686                  type = 2;
687              }
688              int iRow = row[j];
689              if (rhs[iRow] == -COIN_DBL_MAX) {
690                type = 3;
691                break;
692              } else if (rhs[iRow] == -0.5) {
693                // different values
694                unitRhs = false;
695              } else if (rhsValue) {
696                if (rhsValue != rhs[iRow])
697                  unitRhs = false;
698              } else {
699                rhsValue = rhs[iRow];
700              }
701            }
702            // if no elements OK
703            if (type == 3) {
704              // no good
705              possible = 0;
706              break;
707            } else if (type == 2) {
708              if (end - start > 2) {
709                // no good
710                possible = 0;
711                break;
712              } else {
713                // only network
714                possible &= 1;
715                if (!possible)
716                  break;
717              }
718            } else if (type == 1) {
719              // only cover
720              possible &= 2;
721              if (!possible)
722                break;
723            } else if (type == -1) {
724              // only negative cover
725              possible &= 4;
726              if (!possible)
727                break;
728            }
729          }
730        }
731      }
732      if ((possible == 2 || possible == 4) && !unitRhs) {
733#if COIN_DEVELOP > 1
734        printf("XXXXXX Continuous all +1 but different rhs\n");
735#endif
736        possible = 0;
737      }
738      // may be all integer
739      if (possible != 7) {
740        if (!possible)
741          continuousMultiplier = 0.0;
742        else if (possible == 1)
743          continuousMultiplier = 1.0;
744        else
745          continuousMultiplier = 0.0; // 0.5 was incorrect;
746#if COIN_DEVELOP > 1
747        if (continuousMultiplier)
748          printf("XXXXXX multiplier of %g\n", continuousMultiplier);
749#endif
750        if (continuousMultiplier == 0.5) {
751          coeffMultiplier = new double[numberColumns];
752          bool allOne = true;
753          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
754            coeffMultiplier[iColumn] = 1.0;
755            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
756              if (!isInteger(iColumn)) {
757                CoinBigIndex start = columnStart[iColumn];
758                int iRow = row[start];
759                double value = rhs[iRow];
760                assert(value >= 0.0);
761                if (value != 0.0 && value != 1.0)
762                  allOne = false;
763                coeffMultiplier[iColumn] = 0.5 * value;
764              }
765            }
766          }
767          if (allOne) {
768            // back to old way
769            delete[] coeffMultiplier;
770            coeffMultiplier = NULL;
771          }
772        }
773      } else {
774        // all integer
775        problemType_ = problemType;
776#if COIN_DEVELOP > 1
777        printf("Problem type is %d\n", problemType_);
778#endif
779      }
780    }
781
782    // But try again
783    if (continuousMultiplier < 1.0) {
784      memset(rhs, 0, numberRows * sizeof(double));
785      int *count = new int[numberRows];
786      memset(count, 0, numberRows * sizeof(int));
787      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
788        CoinBigIndex start = columnStart[iColumn];
789        CoinBigIndex end = start + columnLength[iColumn];
790        if (upper[iColumn] == lower[iColumn]) {
791          for (CoinBigIndex j = start; j < end; j++) {
792            int iRow = row[j];
793            rhs[iRow] += lower[iColumn] * element[j];
794          }
795        } else if (solver_->isInteger(iColumn)) {
796          for (CoinBigIndex j = start; j < end; j++) {
797            int iRow = row[j];
798            if (fabs(element[j] - floor(element[j] + 0.5)) > 1.0e-10)
799              rhs[iRow] = COIN_DBL_MAX;
800          }
801        } else {
802          for (CoinBigIndex j = start; j < end; j++) {
803            int iRow = row[j];
804            count[iRow]++;
805            if (fabs(element[j]) != 1.0)
806              rhs[iRow] = COIN_DBL_MAX;
807          }
808        }
809      }
810      // now look at continuous
811      bool allGood = true;
812      double direction = solver_->getObjSense();
813      int numberObj = 0;
814      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
815        if (upper[iColumn] > lower[iColumn]) {
816          double objValue = objective[iColumn] * direction;
817          if (objValue && !solver_->isInteger(iColumn)) {
818            numberObj++;
819            CoinBigIndex start = columnStart[iColumn];
820            CoinBigIndex end = start + columnLength[iColumn];
821            if (objValue > 0.0) {
822              // wants to be as low as possible
823              if (lower[iColumn] < -1.0e10 || fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
824                allGood = false;
825                break;
826              } else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
827                allGood = false;
828                break;
829              }
830              bool singletonRow = true;
831              bool equality = false;
832              for (CoinBigIndex j = start; j < end; j++) {
833                int iRow = row[j];
834                if (count[iRow] > 1)
835                  singletonRow = false;
836                else if (rowLower[iRow] == rowUpper[iRow])
837                  equality = true;
838                double rhsValue = rhs[iRow];
839                double lowerValue = rowLower[iRow];
840                double upperValue = rowUpper[iRow];
841                if (rhsValue < 1.0e20) {
842                  if (lowerValue > -1.0e20)
843                    lowerValue -= rhsValue;
844                  if (upperValue < 1.0e20)
845                    upperValue -= rhsValue;
846                }
847                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
848                  || fabs(element[j]) != 1.0) {
849                  // no good
850                  allGood = false;
851                  break;
852                }
853                if (element[j] > 0.0) {
854                  if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
855                    // no good
856                    allGood = false;
857                    break;
858                  }
859                } else {
860                  if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
861                    // no good
862                    allGood = false;
863                    break;
864                  }
865                }
866              }
867              if (!singletonRow && end > start + 1 && !equality)
868                allGood = false;
869              if (!allGood)
870                break;
871            } else {
872              // wants to be as high as possible
873              if (upper[iColumn] > 1.0e10 || fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
874                allGood = false;
875                break;
876              } else if (lower[iColumn] > -1.0e10 && fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
877                allGood = false;
878                break;
879              }
880              bool singletonRow = true;
881              bool equality = false;
882              for (CoinBigIndex j = start; j < end; j++) {
883                int iRow = row[j];
884                if (count[iRow] > 1)
885                  singletonRow = false;
886                else if (rowLower[iRow] == rowUpper[iRow])
887                  equality = true;
888                double rhsValue = rhs[iRow];
889                double lowerValue = rowLower[iRow];
890                double upperValue = rowUpper[iRow];
891                if (rhsValue < 1.0e20) {
892                  if (lowerValue > -1.0e20)
893                    lowerValue -= rhsValue;
894                  if (upperValue < 1.0e20)
895                    upperValue -= rhsValue;
896                }
897                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
898                  || fabs(element[j]) != 1.0) {
899                  // no good
900                  allGood = false;
901                  break;
902                }
903                if (element[j] < 0.0) {
904                  if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
905                    // no good
906                    allGood = false;
907                    break;
908                  }
909                } else {
910                  if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
911                    // no good
912                    allGood = false;
913                    break;
914                  }
915                }
916              }
917              if (!singletonRow && end > start + 1 && !equality)
918                allGood = false;
919              if (!allGood)
920                break;
921            }
922          }
923        }
924      }
925      delete[] count;
926      if (allGood) {
927#if COIN_DEVELOP > 1
928        if (numberObj)
929          printf("YYYY analysis says all continuous with costs will be integer\n");
930#endif
931        continuousMultiplier = 1.0;
932      }
933    }
934    delete[] rhs;
935  }
936  /*
937      Take a first scan to see if there are unfixed continuous variables in the
938      objective.  If so, the minimum objective change could be arbitrarily small.
939      Also pick off the maximum coefficient of an unfixed integer variable.
940
941      If the objective is found to contain only integer variables, set the
942      fathoming discipline to strict.
943    */
944  double maximumCost = 0.0;
945  //double trueIncrement=0.0;
946  int iColumn;
947  int numberColumns = getNumCols();
948  double scaleFactor = 1.0; // due to rhs etc
949  /*
950      Original model did not have integer bounds.
951    */
952  if ((specialOptions_ & 65536) == 0) {
953    /* be on safe side (later look carefully as may be able to
954           to get 0.5 say if bounds are multiples of 0.5 */
955    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
956      if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
957        double value;
958        value = fabs(lower[iColumn]);
959        if (floor(value + 0.5) != value) {
960          scaleFactor = CoinMin(scaleFactor, 0.5);
961          if (floor(2.0 * value + 0.5) != 2.0 * value) {
962            scaleFactor = CoinMin(scaleFactor, 0.25);
963            if (floor(4.0 * value + 0.5) != 4.0 * value) {
964              scaleFactor = 0.0;
965            }
966          }
967        }
968        value = fabs(upper[iColumn]);
969        if (floor(value + 0.5) != value) {
970          scaleFactor = CoinMin(scaleFactor, 0.5);
971          if (floor(2.0 * value + 0.5) != 2.0 * value) {
972            scaleFactor = CoinMin(scaleFactor, 0.25);
973            if (floor(4.0 * value + 0.5) != 4.0 * value) {
974              scaleFactor = 0.0;
975            }
976          }
977        }
978      }
979    }
980  }
981  bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0;
982  if (possibleMultiple) {
983    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
984      if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
985        maximumCost = CoinMax(maximumCost, fabs(objective[iColumn]));
986      }
987    }
988  }
989  setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple);
990  /*
991      If a nontrivial increment is possible, try and figure it out. We're looking
992      for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
993      variables. Since the c<j> might not be integers, try and inflate them
994      sufficiently that they look like integers (and we'll deflate the gcd
995      later).
996
997      2520.0 is used as it is a nice multiple of 2,3,5,7
998    */
999  if (possibleMultiple && maximumCost) {
1000    int increment = 0;
1001    double multiplier = 2520.0;
1002    while (10.0 * multiplier * maximumCost < 1.0e8)
1003      multiplier *= 10.0;
1004    int bigIntegers = 0; // Count of large costs which are integer
1005    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1006      if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
1007        double objValue = fabs(objective[iColumn]);
1008        if (!isInteger(iColumn)) {
1009          if (!coeffMultiplier)
1010            objValue *= continuousMultiplier;
1011          else
1012            objValue *= coeffMultiplier[iColumn];
1013        }
1014        if (objValue) {
1015          double value = objValue * multiplier;
1016          if (value < 2.1e9) {
1017            int nearest = static_cast<int>(floor(value + 0.5));
1018            if (fabs(value - floor(value + 0.5)) > 1.0e-8) {
1019              increment = 0;
1020              break;
1021            } else if (!increment) {
1022              increment = nearest;
1023            } else {
1024              increment = gcd(increment, nearest);
1025            }
1026          } else {
1027            // large value - may still be multiple of 1.0
1028            if (fabs(objValue - floor(objValue + 0.5)) > 1.0e-8) {
1029              increment = 0;
1030              break;
1031            } else {
1032              bigIntegers++;
1033            }
1034          }
1035        }
1036      }
1037    }
1038    delete[] coeffMultiplier;
1039    /*
1040          If the increment beats the current value for objective change, install it.
1041        */
1042    if (increment) {
1043      double value = increment;
1044      double cutoff = getDblParam(CbcModel::CbcCutoffIncrement);
1045      if (bigIntegers) {
1046        // allow for 1.0
1047        increment = gcd(increment, static_cast<int>(multiplier));
1048        value = increment;
1049      }
1050      value /= multiplier;
1051      value *= scaleFactor;
1052      //trueIncrement=CoinMax(cutoff,value);;
1053      if (value * 0.999 > cutoff) {
1054        messageHandler()->message(CBC_INTEGERINCREMENT,
1055          messages())
1056          << value << CoinMessageEol;
1057        setDblParam(CbcModel::CbcCutoffIncrement, CoinMax(value * 0.999, value - 1.0e-4));
1058      }
1059    }
1060  }
1061
1062  return;
1063}
1064
1065/*
1066saveModel called (carved out of) BranchandBound
1067*/
1068void CbcModel::saveModel(OsiSolverInterface *saveSolver, double *checkCutoffForRestart, bool *feasible)
1069{
1070  if (saveSolver && (specialOptions_ & 32768) != 0) {
1071    // See if worth trying reduction
1072    *checkCutoffForRestart = getCutoff();
1073    bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() && (*checkCutoffForRestart < 1.0e20);
1074    int numberColumns = getNumCols();
1075    if (tryNewSearch) {
1076#if CBC_USEFUL_PRINTING > 1
1077      printf("after %d nodes, cutoff %g - looking\n",
1078        numberNodes_, getCutoff());
1079#endif
1080      saveSolver->resolve();
1081      double direction = saveSolver->getObjSense();
1082      double gap = *checkCutoffForRestart - saveSolver->getObjValue() * direction;
1083      double tolerance;
1084      saveSolver->getDblParam(OsiDualTolerance, tolerance);
1085      if (gap <= 0.0)
1086        gap = tolerance;
1087      gap += 100.0 * tolerance;
1088      double integerTolerance = getDblParam(CbcIntegerTolerance);
1089
1090      const double *lower = saveSolver->getColLower();
1091      const double *upper = saveSolver->getColUpper();
1092      const double *solution = saveSolver->getColSolution();
1093      const double *reducedCost = saveSolver->getReducedCost();
1094
1095      int numberFixed = 0;
1096      int numberFixed2 = 0;
1097      for (int i = 0; i < numberIntegers_; i++) {
1098        int iColumn = integerVariable_[i];
1099        double djValue = direction * reducedCost[iColumn];
1100        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1101          if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1102            saveSolver->setColUpper(iColumn, lower[iColumn]);
1103            numberFixed++;
1104          } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1105            saveSolver->setColLower(iColumn, upper[iColumn]);
1106            numberFixed++;
1107          }
1108        } else {
1109          numberFixed2++;
1110        }
1111      }
1112#ifdef COIN_DEVELOP
1113      /*
1114              We're debugging. (specialOptions 1)
1115            */
1116      if ((specialOptions_ & 1) != 0) {
1117        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger();
1118        if (debugger) {
1119          printf("Contains optimal\n");
1120          OsiSolverInterface *temp = saveSolver->clone();
1121          const double *solution = debugger->optimalSolution();
1122          const double *lower = temp->getColLower();
1123          const double *upper = temp->getColUpper();
1124          int n = temp->getNumCols();
1125          for (int i = 0; i < n; i++) {
1126            if (temp->isInteger(i)) {
1127              double value = floor(solution[i] + 0.5);
1128              assert(value >= lower[i] && value <= upper[i]);
1129              temp->setColLower(i, value);
1130              temp->setColUpper(i, value);
1131            }
1132          }
1133          temp->writeMps("reduced_fix");
1134          delete temp;
1135          saveSolver->writeMps("reduced");
1136        } else {
1137          abort();
1138        }
1139      }
1140      printf("Restart could fix %d integers (%d already fixed)\n",
1141        numberFixed + numberFixed2, numberFixed2);
1142#endif
1143      numberFixed += numberFixed2;
1144      if (numberFixed * 20 < numberColumns)
1145        tryNewSearch = false;
1146    }
1147    if (tryNewSearch) {
1148      // back to solver without cuts?
1149      OsiSolverInterface *solver2 = continuousSolver_->clone();
1150      const double *lower = saveSolver->getColLower();
1151      const double *upper = saveSolver->getColUpper();
1152      for (int i = 0; i < numberIntegers_; i++) {
1153        int iColumn = integerVariable_[i];
1154        solver2->setColLower(iColumn, lower[iColumn]);
1155        solver2->setColUpper(iColumn, upper[iColumn]);
1156      }
1157      // swap
1158      delete saveSolver;
1159      saveSolver = solver2;
1160      double *newSolution = new double[numberColumns];
1161      double objectiveValue = *checkCutoffForRestart;
1162      CbcSerendipity heuristic(*this);
1163      if (bestSolution_)
1164        heuristic.setInputSolution(bestSolution_, bestObjective_);
1165      heuristic.setFractionSmall(0.9);
1166      heuristic.setFeasibilityPumpOptions(1008013);
1167      // Use numberNodes to say how many are original rows
1168      heuristic.setNumberNodes(continuousSolver_->getNumRows());
1169#ifdef COIN_DEVELOP
1170      if (continuousSolver_->getNumRows() < saveSolver->getNumRows())
1171        printf("%d rows added ZZZZZ\n",
1172          solver_->getNumRows() - continuousSolver_->getNumRows());
1173#endif
1174      int returnCode = heuristic.smallBranchAndBound(saveSolver,
1175        -1, newSolution,
1176        objectiveValue,
1177        *checkCutoffForRestart, "Reduce");
1178      if (returnCode < 0) {
1179#ifdef COIN_DEVELOP
1180        printf("Restart - not small enough to do search after fixing\n");
1181#endif
1182        delete[] newSolution;
1183      } else {
1184        if ((returnCode & 1) != 0) {
1185          // increment number of solutions so other heuristics can test
1186          numberSolutions_++;
1187          numberHeuristicSolutions_++;
1188          lastHeuristic_ = NULL;
1189          setBestSolution(CBC_ROUNDING, objectiveValue, newSolution);
1190        }
1191        delete[] newSolution;
1192        *feasible = false; // stop search
1193      }
1194#if 0 // probably not needed def CBC_THREAD
1195            if (master_) {
1196                lockThread();
1197                if (parallelMode() > 0) {
1198                    while (master_->waitForThreadsInTree(0)) {
1199                        lockThread();
1200                        double dummyBest;
1201                        tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
1202                        //unlockThread();
1203                    }
1204                }
1205                master_->waitForThreadsInTree(2);
1206                delete master_;
1207                master_ = NULL;
1208                masterThread_ = NULL;
1209            }
1210#endif
1211    }
1212  }
1213}
1214/*
1215Adds integers, called from BranchandBound()
1216*/
1217void CbcModel::AddIntegers()
1218{
1219  int numberColumns = continuousSolver_->getNumCols();
1220  int numberRows = continuousSolver_->getNumRows();
1221  int numberOriginalIntegers = numberIntegers_;
1222  int *del = new int[CoinMax(numberColumns, numberRows)];
1223  int *original = new int[numberColumns];
1224  char *possibleRow = new char[numberRows];
1225  {
1226    const CoinPackedMatrix *rowCopy = continuousSolver_->getMatrixByRow();
1227    const int *column = rowCopy->getIndices();
1228    const int *rowLength = rowCopy->getVectorLengths();
1229    const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
1230    const double *rowLower = continuousSolver_->getRowLower();
1231    const double *rowUpper = continuousSolver_->getRowUpper();
1232    const double *element = rowCopy->getElements();
1233    for (int i = 0; i < numberRows; i++) {
1234      int nLeft = 0;
1235      bool possible = false;
1236      if (rowLower[i] < -1.0e20) {
1237        double value = rowUpper[i];
1238        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1239          possible = true;
1240      } else if (rowUpper[i] > 1.0e20) {
1241        double value = rowLower[i];
1242        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1243          possible = true;
1244      } else {
1245        double value = rowUpper[i];
1246        if (rowLower[i] == rowUpper[i] && fabs(value - floor(value + 0.5)) < 1.0e-8)
1247          possible = true;
1248      }
1249      double allSame = (possible) ? 0.0 : -1.0;
1250      for (CoinBigIndex j = rowStart[i];
1251           j < rowStart[i] + rowLength[i]; j++) {
1252        int iColumn = column[j];
1253        if (continuousSolver_->isInteger(iColumn)) {
1254          if (fabs(element[j]) != 1.0)
1255            possible = false;
1256        } else {
1257          nLeft++;
1258          if (!allSame) {
1259            allSame = fabs(element[j]);
1260          } else if (allSame > 0.0) {
1261            if (allSame != fabs(element[j]))
1262              allSame = -1.0;
1263          }
1264        }
1265      }
1266      if (nLeft == rowLength[i] && allSame > 0.0)
1267        possibleRow[i] = 2;
1268      else if (possible || !nLeft)
1269        possibleRow[i] = 1;
1270      else
1271        possibleRow[i] = 0;
1272    }
1273  }
1274  int nDel = 0;
1275  for (int i = 0; i < numberColumns; i++) {
1276    original[i] = i;
1277    if (continuousSolver_->isInteger(i))
1278      del[nDel++] = i;
1279  }
1280  {
1281    // we must not exclude current best solution (rounding errors)
1282    // also not if large values
1283    const int *row = continuousSolver_->getMatrixByCol()->getIndices();
1284    const CoinBigIndex *columnStart = continuousSolver_->getMatrixByCol()->getVectorStarts();
1285    const int *columnLength = continuousSolver_->getMatrixByCol()->getVectorLengths();
1286    const double *solution = continuousSolver_->getColSolution();
1287    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1288      if (!continuousSolver_->isInteger(iColumn)) {
1289        double value = bestSolution_ ? bestSolution_[iColumn] : 0.0;
1290        double value2 = solution[iColumn];
1291        if (fabs(value - floor(value + 0.5)) > 1.0e-8 || fabs(value2) > 1.0e3) {
1292          CoinBigIndex start = columnStart[iColumn];
1293          CoinBigIndex end = start + columnLength[iColumn];
1294          for (CoinBigIndex j = start; j < end; j++) {
1295            int iRow = row[j];
1296            possibleRow[iRow] = 0;
1297          }
1298        }
1299      }
1300    }
1301  }
1302  int nExtra = 0;
1303  OsiSolverInterface *copy1 = continuousSolver_->clone();
1304  int nPass = 0;
1305  while (nDel && nPass < 10) {
1306    nPass++;
1307    OsiSolverInterface *copy2 = copy1->clone();
1308    int nLeft = 0;
1309    for (int i = 0; i < nDel; i++)
1310      original[del[i]] = -1;
1311    for (int i = 0; i < numberColumns; i++) {
1312      int kOrig = original[i];
1313      if (kOrig >= 0)
1314        original[nLeft++] = kOrig;
1315    }
1316    assert(nLeft == numberColumns - nDel);
1317    copy2->deleteCols(nDel, del);
1318    numberColumns = copy2->getNumCols();
1319    const CoinPackedMatrix *rowCopy = copy2->getMatrixByRow();
1320    numberRows = rowCopy->getNumRows();
1321    const int *column = rowCopy->getIndices();
1322    const int *rowLength = rowCopy->getVectorLengths();
1323    const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
1324    const double *rowLower = copy2->getRowLower();
1325    const double *rowUpper = copy2->getRowUpper();
1326    const double *element = rowCopy->getElements();
1327    const CoinPackedMatrix *columnCopy = copy2->getMatrixByCol();
1328    const int *columnLength = columnCopy->getVectorLengths();
1329    nDel = 0;
1330    // Could do gcd stuff on ones with costs
1331    for (int i = 0; i < numberRows; i++) {
1332      if (!rowLength[i]) {
1333        del[nDel++] = i;
1334      } else if (possibleRow[i]) {
1335        if (rowLength[i] == 1) {
1336          CoinBigIndex k = rowStart[i];
1337          int iColumn = column[k];
1338          if (!copy2->isInteger(iColumn)) {
1339            double mult = 1.0 / fabs(element[k]);
1340            if (rowLower[i] < -1.0e20) {
1341              // treat rhs as multiple of 1 unless elements all same
1342              double value = ((possibleRow[i] == 2) ? rowUpper[i] : 1.0) * mult;
1343              if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1344                del[nDel++] = i;
1345                if (columnLength[iColumn] == 1) {
1346                  copy2->setInteger(iColumn);
1347                  int kOrig = original[iColumn];
1348                  setOptionalInteger(kOrig);
1349                }
1350              }
1351            } else if (rowUpper[i] > 1.0e20) {
1352              // treat rhs as multiple of 1 unless elements all same
1353              double value = ((possibleRow[i] == 2) ? rowLower[i] : 1.0) * mult;
1354              if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1355                del[nDel++] = i;
1356                if (columnLength[iColumn] == 1) {
1357                  copy2->setInteger(iColumn);
1358                  int kOrig = original[iColumn];
1359                  setOptionalInteger(kOrig);
1360                }
1361              }
1362            } else {
1363              // treat rhs as multiple of 1 unless elements all same
1364              double value = ((possibleRow[i] == 2) ? rowUpper[i] : 1.0) * mult;
1365              if (rowLower[i] == rowUpper[i] && fabs(value - floor(value + 0.5)) < 1.0e-8) {
1366                del[nDel++] = i;
1367                copy2->setInteger(iColumn);
1368                int kOrig = original[iColumn];
1369                setOptionalInteger(kOrig);
1370              }
1371            }
1372          }
1373        } else {
1374          // only if all singletons
1375          bool possible = false;
1376          if (rowLower[i] < -1.0e20) {
1377            double value = rowUpper[i];
1378            if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1379              possible = true;
1380          } else if (rowUpper[i] > 1.0e20) {
1381            double value = rowLower[i];
1382            if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1383              possible = true;
1384          } else {
1385            double value = rowUpper[i];
1386            if (rowLower[i] == rowUpper[i] && fabs(value - floor(value + 0.5)) < 1.0e-8)
1387              possible = true;
1388          }
1389          if (possible) {
1390            for (CoinBigIndex j = rowStart[i];
1391                 j < rowStart[i] + rowLength[i]; j++) {
1392              int iColumn = column[j];
1393              if (columnLength[iColumn] != 1 || fabs(element[j]) != 1.0) {
1394                possible = false;
1395                break;
1396              }
1397            }
1398            if (possible) {
1399              for (CoinBigIndex j = rowStart[i];
1400                   j < rowStart[i] + rowLength[i]; j++) {
1401                int iColumn = column[j];
1402                if (!copy2->isInteger(iColumn)) {
1403                  copy2->setInteger(iColumn);
1404                  int kOrig = original[iColumn];
1405                  setOptionalInteger(kOrig);
1406                }
1407              }
1408              del[nDel++] = i;
1409            }
1410          }
1411        }
1412      }
1413    }
1414    if (nDel) {
1415      copy2->deleteRows(nDel, del);
1416      // pack down possible
1417      int n = 0;
1418      for (int i = 0; i < nDel; i++)
1419        possibleRow[del[i]] = -1;
1420      for (int i = 0; i < numberRows; i++) {
1421        if (possibleRow[i] >= 0)
1422          possibleRow[n++] = possibleRow[i];
1423      }
1424    }
1425    if (nDel != numberRows) {
1426      nDel = 0;
1427      for (int i = 0; i < numberColumns; i++) {
1428        if (copy2->isInteger(i)) {
1429          del[nDel++] = i;
1430          nExtra++;
1431        }
1432      }
1433    } else {
1434      nDel = 0;
1435    }
1436    delete copy1;
1437    copy1 = copy2->clone();
1438    delete copy2;
1439  }
1440  // See if what's left is a network
1441  bool couldBeNetwork = false;
1442  if (copy1->getNumRows() && copy1->getNumCols()) {
1443#ifdef COIN_HAS_CLP
1444    OsiClpSolverInterface *clpSolver
1445      = dynamic_cast<OsiClpSolverInterface *>(copy1);
1446    if (false && clpSolver) {
1447      numberRows = clpSolver->getNumRows();
1448      char *rotate = new char[numberRows];
1449      int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0);
1450      delete[] rotate;
1451#if CBC_USEFUL_PRINTING > 1
1452      printf("INTA network %d rows out of %d\n", n, numberRows);
1453#endif
1454      if (CoinAbs(n) == numberRows) {
1455        couldBeNetwork = true;
1456        for (int i = 0; i < numberRows; i++) {
1457          if (!possibleRow[i]) {
1458            couldBeNetwork = false;
1459#if CBC_USEFUL_PRINTING > 1
1460            printf("but row %d is bad\n", i);
1461#endif
1462            break;
1463          }
1464        }
1465      }
1466    } else
1467#endif
1468    {
1469      numberColumns = copy1->getNumCols();
1470      numberRows = copy1->getNumRows();
1471      const double *rowLower = copy1->getRowLower();
1472      const double *rowUpper = copy1->getRowUpper();
1473      couldBeNetwork = true;
1474      for (int i = 0; i < numberRows; i++) {
1475        if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) {
1476          couldBeNetwork = false;
1477          break;
1478        }
1479        if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) {
1480          couldBeNetwork = false;
1481          break;
1482        }
1483        if (possibleRow[i] == 0) {
1484          couldBeNetwork = false;
1485          break;
1486        }
1487      }
1488      if (couldBeNetwork) {
1489        const CoinPackedMatrix *matrixByCol = copy1->getMatrixByCol();
1490        const double *element = matrixByCol->getElements();
1491        //const int * row = matrixByCol->getIndices();
1492        const CoinBigIndex *columnStart = matrixByCol->getVectorStarts();
1493        const int *columnLength = matrixByCol->getVectorLengths();
1494        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1495          CoinBigIndex start = columnStart[iColumn];
1496          CoinBigIndex end = start + columnLength[iColumn];
1497          if (end > start + 2) {
1498            couldBeNetwork = false;
1499            break;
1500          }
1501          int type = 0;
1502          for (CoinBigIndex j = start; j < end; j++) {
1503            double value = element[j];
1504            if (fabs(value) != 1.0) {
1505              couldBeNetwork = false;
1506              break;
1507            } else if (value == 1.0) {
1508              if ((type & 1) == 0)
1509                type |= 1;
1510              else
1511                type = 7;
1512            } else if (value == -1.0) {
1513              if ((type & 2) == 0)
1514                type |= 2;
1515              else
1516                type = 7;
1517            }
1518          }
1519          if (type > 3) {
1520            couldBeNetwork = false;
1521            break;
1522          }
1523        }
1524      }
1525    }
1526  }
1527  if (couldBeNetwork) {
1528    for (int i = 0; i < numberColumns; i++)
1529      setOptionalInteger(original[i]);
1530  }
1531  if (nExtra || couldBeNetwork) {
1532    numberColumns = copy1->getNumCols();
1533    numberRows = copy1->getNumRows();
1534    if (!numberColumns || !numberRows) {
1535      int numberColumns = solver_->getNumCols();
1536      for (int i = 0; i < numberColumns; i++)
1537        assert(solver_->isInteger(i));
1538    }
1539#if CBC_USEFUL_PRINTING > 1
1540    if (couldBeNetwork || nExtra)
1541      printf("INTA %d extra integers, %d left%s\n", nExtra,
1542        numberColumns,
1543        couldBeNetwork ? ", all network" : "");
1544#endif
1545    findIntegers(true, 2);
1546    convertToDynamic();
1547  }
1548#if CBC_USEFUL_PRINTING > 1
1549  if (!couldBeNetwork && copy1->getNumCols() && copy1->getNumRows()) {
1550    printf("INTA %d rows and %d columns remain\n",
1551      copy1->getNumRows(), copy1->getNumCols());
1552    if (copy1->getNumCols() < 200) {
1553      copy1->writeMps("moreint");
1554      printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n",
1555        copy1->getNumRows(), copy1->getNumCols());
1556    }
1557  }
1558#endif
1559  delete copy1;
1560  delete[] del;
1561  delete[] original;
1562  delete[] possibleRow;
1563  // double check increment
1564  analyzeObjective();
1565  // If any changes - tell code
1566  if (numberOriginalIntegers < numberIntegers_)
1567    synchronizeModel();
1568}
1569/**
1570  \todo
1571  Normally, it looks like we enter here from command dispatch in the main
1572  routine, after calling the solver for an initial solution
1573  (CbcModel::initialSolve, which simply calls the solver's initialSolve
1574  routine.) The first thing we do is call resolve. Presumably there are
1575  circumstances where this is nontrivial? There's also a call from
1576  CbcModel::originalModel (tied up with integer presolve), which should be
1577  checked.
1578
1579*/
1580
1581#ifdef CONFLICT_CUTS
1582#if PRINT_CONFLICT == 1
1583static int numberConflictCuts = 0;
1584static int lastNumberConflictCuts = 0;
1585static double lengthConflictCuts = 0.0;
1586#endif
1587#endif
1588/*
1589  The overall flow can be divided into three stages:
1590    * Prep: Check that the lp relaxation remains feasible at the root. If so,
1591      do all the setup for B&C.
1592    * Process the root node: Generate cuts, apply heuristics, and in general do
1593      the best we can to resolve the problem without B&C.
1594    * Do B&C search until we hit a limit or exhaust the search tree.
1595
1596  Keep in mind that in general there is no node in the search tree that
1597  corresponds to the active subproblem. The active subproblem is represented
1598  by the current state of the model,  of the solver, and of the constraint
1599  system held by the solver.
1600*/
1601#ifdef COIN_HAS_CPX
1602#include "OsiCpxSolverInterface.hpp"
1603#include "cplex.h"
1604#endif
1605void CbcModel::branchAndBound(int doStatistics)
1606
1607{
1608  if (!parentModel_) {
1609    /*
1610        Capture a time stamp before we start (unless set).
1611      */
1612    if (!dblParam_[CbcStartSeconds]) {
1613      if (!useElapsedTime())
1614        dblParam_[CbcStartSeconds] = CoinCpuTime();
1615      else
1616        dblParam_[CbcStartSeconds] = CoinGetTimeOfDay();
1617    }
1618  }
1619  dblParam_[CbcSmallestChange] = COIN_DBL_MAX;
1620  dblParam_[CbcSumChange] = 0.0;
1621  dblParam_[CbcLargestChange] = 0.0;
1622  intParam_[CbcNumberBranches] = 0;
1623  double lastBestPossibleObjective = -COIN_DBL_MAX;
1624  // when to check for restart
1625  int nextCheckRestart = 50;
1626  // Force minimization !!!!
1627  bool flipObjective = (solver_->getObjSense() < 0.0);
1628  if (flipObjective)
1629    flipModel();
1630  dblParam_[CbcOptimizationDirection] = 1.0; // was solver_->getObjSense();
1631  strongInfo_[0] = 0;
1632  strongInfo_[1] = 0;
1633  strongInfo_[2] = 0;
1634  strongInfo_[3] = 0;
1635  strongInfo_[4] = 0;
1636  strongInfo_[5] = 0;
1637  strongInfo_[6] = 0;
1638  numberStrongIterations_ = 0;
1639  currentNode_ = NULL;
1640  // See if should do cuts old way
1641  if (parallelMode() < 0) {
1642    specialOptions_ |= 4096 + 8192;
1643  } else if (parallelMode() > 0) {
1644    specialOptions_ |= 4096;
1645  }
1646  int saveMoreSpecialOptions = moreSpecialOptions_;
1647  if (dynamic_cast<CbcTreeLocal *>(tree_))
1648    specialOptions_ |= 4096 + 8192;
1649#ifdef COIN_HAS_CLP
1650  {
1651    OsiClpSolverInterface *clpSolver
1652      = dynamic_cast<OsiClpSolverInterface *>(solver_);
1653    if (clpSolver) {
1654      // pass in disaster handler
1655      CbcDisasterHandler handler(this);
1656      clpSolver->passInDisasterHandler(&handler);
1657      // Initialise solvers seed (unless users says not)
1658      if ((specialOptions_ & 4194304) == 0)
1659        clpSolver->getModelPtr()->setRandomSeed(1234567);
1660#ifdef JJF_ZERO
1661      // reduce factorization frequency
1662      int frequency = clpSolver->getModelPtr()->factorizationFrequency();
1663      clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
1664#endif
1665    }
1666  }
1667#endif
1668  // original solver (only set if pre-processing)
1669  OsiSolverInterface *originalSolver = NULL;
1670  int numberOriginalObjects = numberObjects_;
1671  OsiObject **originalObject = NULL;
1672  // Save whether there were any objects
1673  bool noObjects = (numberObjects_ == 0);
1674  // Set up strategies
1675  /*
1676      See if the user has supplied a strategy object and deal with it if present.
1677      The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
1678      perform integer preprocessing, if requested.
1679
1680      We need to hang on to a pointer to solver_. setupOther will assign a
1681      preprocessed solver to model, but will instruct assignSolver not to trash the
1682      existing one.
1683    */
1684  if (strategy_) {
1685    // May do preprocessing
1686    originalSolver = solver_;
1687    strategy_->setupOther(*this);
1688    if (strategy_->preProcessState()) {
1689      // pre-processing done
1690      if (strategy_->preProcessState() < 0) {
1691        // infeasible (or unbounded)
1692        status_ = 0;
1693        if (!solver_->isProvenDualInfeasible()) {
1694          handler_->message(CBC_INFEAS, messages_) << CoinMessageEol;
1695          secondaryStatus_ = 1;
1696        } else {
1697          handler_->message(CBC_UNBOUNDED,
1698            messages_)
1699            << CoinMessageEol;
1700          secondaryStatus_ = 7;
1701        }
1702        originalContinuousObjective_ = COIN_DBL_MAX;
1703        if (flipObjective)
1704          flipModel();
1705        return;
1706      } else if (numberObjects_ && object_) {
1707        numberOriginalObjects = numberObjects_;
1708        // redo sequence
1709        numberIntegers_ = 0;
1710        int numberColumns = getNumCols();
1711        int nOrig = originalSolver->getNumCols();
1712        CglPreProcess *process = strategy_->process();
1713        assert(process);
1714        const int *originalColumns = process->originalColumns();
1715        // allow for cliques etc
1716        nOrig = CoinMax(nOrig, originalColumns[numberColumns - 1] + 1);
1717        // try and redo debugger
1718        OsiRowCutDebugger *debugger = const_cast<OsiRowCutDebugger *>(solver_->getRowCutDebuggerAlways());
1719        if (debugger) {
1720          if (numberColumns <= debugger->numberColumns())
1721            debugger->redoSolution(numberColumns, originalColumns);
1722          else
1723            debugger = NULL; // no idea how to handle (SOS?)
1724        }
1725        // User-provided solution might have been best. Synchronise.
1726        if (bestSolution_) {
1727          // need to redo - in case no better found in BAB
1728          // just get integer part right
1729          for (int i = 0; i < numberColumns; i++) {
1730            int jColumn = originalColumns[i];
1731            bestSolution_[i] = bestSolution_[jColumn];
1732          }
1733        }
1734        originalObject = object_;
1735        // object number or -1
1736        int *temp = new int[nOrig];
1737        int iColumn;
1738        for (iColumn = 0; iColumn < nOrig; iColumn++)
1739          temp[iColumn] = -1;
1740        int iObject;
1741        int nNonInt = 0;
1742        for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1743          iColumn = originalObject[iObject]->columnNumber();
1744          if (iColumn < 0) {
1745            nNonInt++;
1746          } else {
1747            temp[iColumn] = iObject;
1748          }
1749        }
1750        int numberNewIntegers = 0;
1751        int numberOldIntegers = 0;
1752        int numberOldOther = 0;
1753        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1754          int jColumn = originalColumns[iColumn];
1755          if (temp[jColumn] >= 0) {
1756            int iObject = temp[jColumn];
1757            CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(originalObject[iObject]);
1758            if (obj)
1759              numberOldIntegers++;
1760            else
1761              numberOldOther++;
1762          } else if (isInteger(iColumn)) {
1763            numberNewIntegers++;
1764          }
1765        }
1766        /*
1767                  Allocate an array to hold the indices of the integer variables.
1768                  Make a large enough array for all objects
1769                */
1770        numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
1771        object_ = new OsiObject *[numberObjects_];
1772        delete[] integerVariable_;
1773        integerVariable_ = new int[numberNewIntegers + numberOldIntegers];
1774        /*
1775                  Walk the variables again, filling in the indices and creating objects for
1776                  the integer variables. Initially, the objects hold the index and upper &
1777                  lower bounds.
1778                */
1779        numberIntegers_ = 0;
1780        int n = originalColumns[numberColumns - 1] + 1;
1781        int *backward = new int[n];
1782        int i;
1783        for (i = 0; i < n; i++)
1784          backward[i] = -1;
1785        for (i = 0; i < numberColumns; i++)
1786          backward[originalColumns[i]] = i;
1787        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1788          int jColumn = originalColumns[iColumn];
1789          if (temp[jColumn] >= 0) {
1790            int iObject = temp[jColumn];
1791            CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(originalObject[iObject]);
1792            if (obj) {
1793              object_[numberIntegers_] = originalObject[iObject]->clone();
1794              // redo ids etc
1795              //object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
1796              object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
1797              integerVariable_[numberIntegers_++] = iColumn;
1798            }
1799          } else if (isInteger(iColumn)) {
1800            object_[numberIntegers_] = new CbcSimpleInteger(this, iColumn);
1801            integerVariable_[numberIntegers_++] = iColumn;
1802          }
1803        }
1804        delete[] backward;
1805        numberObjects_ = numberIntegers_;
1806        // Now append other column stuff
1807        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1808          int jColumn = originalColumns[iColumn];
1809          if (temp[jColumn] >= 0) {
1810            int iObject = temp[jColumn];
1811            CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(originalObject[iObject]);
1812            if (!obj) {
1813              object_[numberObjects_] = originalObject[iObject]->clone();
1814              // redo ids etc
1815              CbcObject *obj = dynamic_cast<CbcObject *>(object_[numberObjects_]);
1816              assert(obj);
1817              obj->redoSequenceEtc(this, numberColumns, originalColumns);
1818              numberObjects_++;
1819            }
1820          }
1821        }
1822        // now append non column stuff
1823        for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1824          iColumn = originalObject[iObject]->columnNumber();
1825          if (iColumn < 0) {
1826            // already has column numbers changed
1827            object_[numberObjects_] = originalObject[iObject]->clone();
1828#ifdef JJF_ZERO
1829            // redo ids etc
1830            CbcObject *obj = dynamic_cast<CbcObject *>(object_[numberObjects_]);
1831            assert(obj);
1832            obj->redoSequenceEtc(this, numberColumns, originalColumns);
1833#endif
1834            numberObjects_++;
1835          }
1836        }
1837        delete[] temp;
1838        if (!numberObjects_)
1839          handler_->message(CBC_NOINT, messages_) << CoinMessageEol;
1840      } else {
1841        int numberColumns = getNumCols();
1842        CglPreProcess *process = strategy_->process();
1843        assert(process);
1844        const int *originalColumns = process->originalColumns();
1845        // try and redo debugger
1846        OsiRowCutDebugger *debugger = const_cast<OsiRowCutDebugger *>(solver_->getRowCutDebuggerAlways());
1847        if (debugger)
1848          debugger->redoSolution(numberColumns, originalColumns);
1849      }
1850    } else {
1851      //no preprocessing
1852      originalSolver = NULL;
1853    }
1854    strategy_->setupCutGenerators(*this);
1855    strategy_->setupHeuristics(*this);
1856    // Set strategy print level to models
1857    strategy_->setupPrinting(*this, handler_->logLevel());
1858  }
1859  eventHappened_ = false;
1860  CbcEventHandler *eventHandler = getEventHandler();
1861  if (eventHandler)
1862    eventHandler->setModel(this);
1863#define CLIQUE_ANALYSIS
1864#ifdef CLIQUE_ANALYSIS
1865  // set up for probing
1866  // If we're doing clever stuff with cliques, additional info here.
1867  if (!parentModel_)
1868    probingInfo_ = new CglTreeProbingInfo(solver_);
1869  else
1870    probingInfo_ = NULL;
1871#else
1872  probingInfo_ = NULL;
1873#endif
1874
1875  // Try for dominated columns
1876  if ((specialOptions_ & 64) != 0) {
1877    CglDuplicateRow dupcuts(solver_);
1878    dupcuts.setMode(2);
1879    CglStored *storedCuts = dupcuts.outDuplicates(solver_);
1880    if (storedCuts) {
1881      COIN_DETAIL_PRINT(printf("adding dup cuts\n"));
1882      addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
1883        true, false, false, -200);
1884    }
1885  }
1886  if (!nodeCompare_)
1887    nodeCompare_ = new CbcCompareDefault();
1888  ;
1889  // See if hot start wanted
1890  CbcCompareBase *saveCompare = NULL;
1891  // User supplied hotstart. Adapt for preprocessing.
1892  if (hotstartSolution_) {
1893    if (strategy_ && strategy_->preProcessState() > 0) {
1894      CglPreProcess *process = strategy_->process();
1895      assert(process);
1896      int n = solver_->getNumCols();
1897      const int *originalColumns = process->originalColumns();
1898      // columns should be in order ... but
1899      double *tempS = new double[n];
1900      for (int i = 0; i < n; i++) {
1901        int iColumn = originalColumns[i];
1902        tempS[i] = hotstartSolution_[iColumn];
1903      }
1904      delete[] hotstartSolution_;
1905      hotstartSolution_ = tempS;
1906      if (hotstartPriorities_) {
1907        int *tempP = new int[n];
1908        for (int i = 0; i < n; i++) {
1909          int iColumn = originalColumns[i];
1910          tempP[i] = hotstartPriorities_[iColumn];
1911        }
1912        delete[] hotstartPriorities_;
1913        hotstartPriorities_ = tempP;
1914      }
1915    }
1916    saveCompare = nodeCompare_;
1917    // depth first
1918    nodeCompare_ = new CbcCompareDepth();
1919  }
1920  if (!problemFeasibility_)
1921    problemFeasibility_ = new CbcFeasibilityBase();
1922#ifdef CBC_DEBUG
1923  std::string problemName;
1924  solver_->getStrParam(OsiProbName, problemName);
1925  printf("Problem name - %s\n", problemName.c_str());
1926  solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0);
1927#endif
1928  /*
1929      Assume we're done, and see if we're proven wrong.
1930    */
1931  status_ = 0;
1932  secondaryStatus_ = 0;
1933  phase_ = 0;
1934  /*
1935      Scan the variables, noting the integer variables. Create an
1936      CbcSimpleInteger object for each integer variable.
1937    */
1938  findIntegers(false);
1939  // Say not dynamic pseudo costs
1940  ownership_ &= ~0x40000000;
1941  // If dynamic pseudo costs then do
1942  if (numberBeforeTrust_)
1943    convertToDynamic();
1944  // Set up char array to say if integer (speed)
1945  delete[] integerInfo_;
1946  {
1947    int n = solver_->getNumCols();
1948    integerInfo_ = new char[n];
1949    for (int i = 0; i < n; i++) {
1950      if (solver_->isInteger(i))
1951        integerInfo_[i] = 1;
1952      else
1953        integerInfo_[i] = 0;
1954    }
1955  }
1956  if (preferredWay_) {
1957    // set all unset ones
1958    for (int iObject = 0; iObject < numberObjects_; iObject++) {
1959      CbcObject *obj = dynamic_cast<CbcObject *>(object_[iObject]);
1960      if (obj && !obj->preferredWay())
1961        obj->setPreferredWay(preferredWay_);
1962    }
1963  }
1964  /*
1965      Ensure that objects on the lists of OsiObjects, heuristics, and cut
1966      generators attached to this model all refer to this model.
1967    */
1968  synchronizeModel();
1969  if (!solverCharacteristics_) {
1970    OsiBabSolver *solverCharacteristics = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
1971    if (solverCharacteristics) {
1972      solverCharacteristics_ = solverCharacteristics;
1973    } else {
1974      // replace in solver
1975      OsiBabSolver defaultC;
1976      solver_->setAuxiliaryInfo(&defaultC);
1977      solverCharacteristics_ = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
1978    }
1979  }
1980
1981  solverCharacteristics_->setSolver(solver_);
1982  // Set so we can tell we are in initial phase in resolve
1983  continuousObjective_ = -COIN_DBL_MAX;
1984  /*
1985      Solve the relaxation.
1986
1987      Apparently there are circumstances where this will be non-trivial --- i.e.,
1988      we've done something since initialSolve that's trashed the solution to the
1989      continuous relaxation.
1990    */
1991  /* Tell solver we are in Branch and Cut
1992       Could use last parameter for subtle differences */
1993  solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL);
1994#ifdef COIN_HAS_CLP
1995  {
1996    OsiClpSolverInterface *clpSolver
1997      = dynamic_cast<OsiClpSolverInterface *>(solver_);
1998    if (clpSolver) {
1999      ClpSimplex *clpSimplex = clpSolver->getModelPtr();
2000      if ((specialOptions_ & 32) == 0) {
2001        // take off names (unless going to be saving)
2002        if (numberAnalyzeIterations_ >= 0 || (-numberAnalyzeIterations_ & 64) == 0)
2003          clpSimplex->dropNames();
2004      }
2005      // no crunch if mostly continuous
2006      if ((clpSolver->specialOptions() & (1 + 8)) != (1 + 8)) {
2007        int numberColumns = solver_->getNumCols();
2008        if (numberColumns > 1000 && numberIntegers_ * 4 < numberColumns)
2009          clpSolver->setSpecialOptions(clpSolver->specialOptions() & (~1));
2010      }
2011      //#define NO_CRUNCH
2012#ifdef NO_CRUNCH
2013      printf("TEMP switching off crunch\n");
2014      int iOpt = clpSolver->specialOptions();
2015      iOpt &= ~1;
2016      iOpt |= 65536;
2017      clpSolver->setSpecialOptions(iOpt);
2018#endif
2019    }
2020  }
2021#endif
2022  bool feasible;
2023  numberSolves_ = 0;
2024  {
2025    // check
2026    int numberOdd = 0;
2027    for (int i = 0; i < numberObjects_; i++) {
2028      CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(object_[i]);
2029      if (!obj)
2030        numberOdd++;
2031    }
2032    if (numberOdd)
2033      moreSpecialOptions_ |= 1073741824;
2034  }
2035  // If NLP then we assume already solved outside branchAndbound
2036  if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
2037    feasible = resolve(NULL, 0) != 0;
2038  } else {
2039    // pick up given status
2040    feasible = (solver_->isProvenOptimal() && !solver_->isDualObjectiveLimitReached());
2041  }
2042  if (problemFeasibility_->feasible(this, 0) < 0) {
2043    feasible = false; // pretend infeasible
2044  }
2045  numberSavedSolutions_ = 0;
2046  int saveNumberStrong = numberStrong_;
2047  int saveNumberBeforeTrust = numberBeforeTrust_;
2048  /*
2049      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2050      continue with processing the root node.
2051    */
2052  if (!feasible) {
2053    status_ = 0;
2054    if (!solver_->isProvenDualInfeasible()) {
2055      handler_->message(CBC_INFEAS, messages_) << CoinMessageEol;
2056      secondaryStatus_ = 1;
2057    } else {
2058      handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol;
2059      secondaryStatus_ = 7;
2060    }
2061    originalContinuousObjective_ = COIN_DBL_MAX;
2062    if (bestSolution_ && ((specialOptions_ & 8388608) == 0 || (specialOptions_ & 2048) != 0)) {
2063      // best solution found by various heuristics - set solution
2064      char general[200];
2065      sprintf(general, "Solution of %g already found by heuristic",
2066        bestObjective_);
2067      messageHandler()->message(CBC_GENERAL,
2068        messages())
2069        << general << CoinMessageEol;
2070      setCutoff(1.0e50); // As best solution should be worse than cutoff
2071      // change cutoff as constraint if wanted
2072      if (cutoffRowNumber_ >= 0) {
2073        if (solver_->getNumRows() > cutoffRowNumber_)
2074          solver_->setRowUpper(cutoffRowNumber_, 1.0e50);
2075      }
2076      // also in continuousSolver_
2077      if (continuousSolver_) {
2078        // Solvers know about direction
2079        double direction = solver_->getObjSense();
2080        continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
2081      } else {
2082        continuousSolver_ = solver_->clone();
2083      }
2084      phase_ = 5;
2085      double increment = getDblParam(CbcModel::CbcCutoffIncrement);
2086      if ((specialOptions_ & 4) == 0)
2087        bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
2088      setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1);
2089      continuousSolver_->resolve();
2090      if (!continuousSolver_->isProvenOptimal()) {
2091        continuousSolver_->messageHandler()->setLogLevel(2);
2092        continuousSolver_->initialSolve();
2093      }
2094      delete solver_;
2095      solverCharacteristics_ = NULL;
2096      solver_ = continuousSolver_;
2097      setPointers(solver_);
2098      continuousSolver_ = NULL;
2099    }
2100    solverCharacteristics_ = NULL;
2101    if (flipObjective)
2102      flipModel();
2103    return;
2104  } else if (!numberObjects_) {
2105    // nothing to do
2106    // Undo preprocessing performed during BaB.
2107    if (strategy_ && strategy_->preProcessState() > 0) {
2108      // undo preprocessing
2109      CglPreProcess *process = strategy_->process();
2110      assert(process);
2111      int n = originalSolver->getNumCols();
2112      if (bestSolution_) {
2113        delete[] bestSolution_;
2114        bestSolution_ = new double[n];
2115        process->postProcess(*solver_);
2116      }
2117      strategy_->deletePreProcess();
2118      // Solution now back in originalSolver
2119      delete solver_;
2120      solver_ = originalSolver;
2121      if (bestSolution_) {
2122        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2123        memcpy(bestSolution_, solver_->getColSolution(), n * sizeof(double));
2124      }
2125      // put back original objects if there were any
2126      if (originalObject) {
2127        int iColumn;
2128        assert(ownObjects_);
2129        for (iColumn = 0; iColumn < numberObjects_; iColumn++)
2130          delete object_[iColumn];
2131        delete[] object_;
2132        numberObjects_ = numberOriginalObjects;
2133        object_ = originalObject;
2134        delete[] integerVariable_;
2135        numberIntegers_ = 0;
2136        for (iColumn = 0; iColumn < n; iColumn++) {
2137          if (solver_->isInteger(iColumn))
2138            numberIntegers_++;
2139        }
2140        integerVariable_ = new int[numberIntegers_];
2141        numberIntegers_ = 0;
2142        for (iColumn = 0; iColumn < n; iColumn++) {
2143          if (solver_->isInteger(iColumn))
2144            integerVariable_[numberIntegers_++] = iColumn;
2145        }
2146      }
2147    }
2148    if (flipObjective)
2149      flipModel();
2150    solverCharacteristics_ = NULL;
2151    bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2152    int numberColumns = solver_->getNumCols();
2153    delete[] bestSolution_;
2154    bestSolution_ = new double[numberColumns];
2155    CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2156    return;
2157  }
2158  /*
2159      See if we're using the Osi side of the branching hierarchy. If so, either
2160      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2161      first case, CbcModel owns the objects on the object_ list. In the second
2162      case, the solver holds the objects and object_ simply points to the
2163      solver's list.
2164
2165      080417 The conversion code here (the block protected by `if (obj)') cannot
2166      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2167      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2168      -> all other Cbc object classes. It's structurally impossible for any Osi
2169      object to descend from CbcObject. The only thing I can see is that this is
2170      really dead code, and object detection is now handled from the Osi side.
2171    */
2172  // Convert to Osi if wanted
2173  //OsiBranchingInformation * persistentInfo = NULL;
2174  if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2175    //persistentInfo = new OsiBranchingInformation(solver_);
2176    if (numberOriginalObjects) {
2177      for (int iObject = 0; iObject < numberObjects_; iObject++) {
2178        CbcObject *obj = dynamic_cast<CbcObject *>(object_[iObject]);
2179        if (obj) {
2180          CbcSimpleInteger *obj2 = dynamic_cast<CbcSimpleInteger *>(obj);
2181          if (obj2) {
2182            // back to Osi land
2183            object_[iObject] = obj2->osiObject();
2184            delete obj;
2185          } else {
2186            OsiSimpleInteger *obj3 = dynamic_cast<OsiSimpleInteger *>(obj);
2187            if (!obj3) {
2188              OsiSOS *obj4 = dynamic_cast<OsiSOS *>(obj);
2189              if (!obj4) {
2190                CbcSOS *obj5 = dynamic_cast<CbcSOS *>(obj);
2191                if (obj5) {
2192                  // back to Osi land
2193                  object_[iObject] = obj5->osiObject(solver_);
2194                } else {
2195                  printf("Code up CbcObject type in Osi land\n");
2196                  abort();
2197                }
2198              }
2199            }
2200          }
2201        }
2202      }
2203      // and add to solver
2204      //if (!solver_->numberObjects()) {
2205      solver_->addObjects(numberObjects_, object_);
2206      //} else {
2207      //if (solver_->numberObjects()!=numberOriginalObjects) {
2208      //printf("should have trapped that solver has objects before\n");
2209      //abort();
2210      //}
2211      //}
2212    } else {
2213      /*
2214              As of 080104, findIntegersAndSOS is misleading --- the default OSI
2215              implementation finds only integers.
2216            */
2217      // do from solver
2218      deleteObjects(false);
2219      solver_->findIntegersAndSOS(false);
2220      numberObjects_ = solver_->numberObjects();
2221      object_ = solver_->objects();
2222      ownObjects_ = false;
2223    }
2224    branchingMethod_->chooseMethod()->setSolver(solver_);
2225  }
2226  // take off heuristics if have to (some do not work with SOS, for example)
2227  // object should know what's safe.
2228  {
2229    int numberOdd = 0;
2230    int numberSOS = 0;
2231    for (int i = 0; i < numberObjects_; i++) {
2232      if (!object_[i]->canDoHeuristics())
2233        numberOdd++;
2234      CbcSOS *obj = dynamic_cast<CbcSOS *>(object_[i]);
2235      if (obj)
2236        numberSOS++;
2237    }
2238    if (numberOdd) {
2239      if (numberHeuristics_ && (specialOptions_ & 1024) == 0) {
2240        int k = 0;
2241        for (int i = 0; i < numberHeuristics_; i++) {
2242          if (!heuristic_[i]->canDealWithOdd())
2243            delete heuristic_[i];
2244          else
2245            heuristic_[k++] = heuristic_[i];
2246        }
2247        if (!k) {
2248          delete[] heuristic_;
2249          heuristic_ = NULL;
2250        }
2251        numberHeuristics_ = k;
2252        handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol;
2253      }
2254      // If odd switch off AddIntegers
2255      specialOptions_ &= ~65536;
2256      // switch off fast nodes for now
2257      fastNodeDepth_ = -1;
2258      moreSpecialOptions_ &= ~33554432; // no diving
2259    } else if (numberSOS) {
2260      specialOptions_ |= 128; // say can do SOS in dynamic mode
2261      // switch off fast nodes for now
2262      fastNodeDepth_ = -1;
2263      moreSpecialOptions_ &= ~33554432; // no diving
2264    }
2265    if (numberThreads_ > 0) {
2266      /* switch off fast nodes for now
2267               Trouble is that by time mini bab finishes code is
2268               looking at a different node
2269             */
2270      fastNodeDepth_ = -1;
2271    }
2272  }
2273  // Save objective (just so user can access it)
2274  originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
2275  bestPossibleObjective_ = originalContinuousObjective_;
2276  sumChangeObjective1_ = 0.0;
2277  sumChangeObjective2_ = 0.0;
2278  /*
2279      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2280      Assuming it recognises the problem, when called upon it will check a cut to
2281      see if it cuts off the optimal answer.
2282    */
2283  // If debugger exists set specialOptions_ bit
2284  if (solver_->getRowCutDebuggerAlways()) {
2285    specialOptions_ |= 1;
2286  }
2287
2288#ifdef CBC_DEBUG
2289  if ((specialOptions_ & 1) == 0)
2290    solver_->activateRowCutDebugger(problemName.c_str());
2291  if (solver_->getRowCutDebuggerAlways())
2292    specialOptions_ |= 1;
2293#endif
2294
2295  /*
2296      Begin setup to process a feasible root node.
2297    */
2298  bestObjective_ = CoinMin(bestObjective_, 1.0e50);
2299  if (!bestSolution_) {
2300    numberSolutions_ = 0;
2301    numberHeuristicSolutions_ = 0;
2302  }
2303  stateOfSearch_ = 0;
2304  // Everything is minimization
2305  {
2306    // needed to sync cutoffs
2307    double value;
2308    solver_->getDblParam(OsiDualObjectiveLimit, value);
2309    dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2310  }
2311  double cutoff = getCutoff();
2312  double direction = solver_->getObjSense();
2313  dblParam_[CbcOptimizationDirection] = direction;
2314  if (cutoff < 1.0e20 && direction < 0.0)
2315    messageHandler()->message(CBC_CUTOFF_WARNING1,
2316      messages())
2317      << cutoff << -cutoff << CoinMessageEol;
2318  if (cutoff > bestObjective_)
2319    cutoff = bestObjective_;
2320  setCutoff(cutoff);
2321  /*
2322      We probably already have a current solution, but just in case ...
2323    */
2324  int numberColumns = getNumCols();
2325  if (!currentSolution_)
2326    currentSolution_ = new double[numberColumns];
2327  testSolution_ = currentSolution_;
2328  /*
2329      Create a copy of the solver, thus capturing the original (root node)
2330      constraint system (aka the continuous system).
2331    */
2332  delete continuousSolver_;
2333  continuousSolver_ = solver_->clone();
2334#ifdef CONFLICT_CUTS
2335  if ((moreSpecialOptions_ & 4194304) != 0) {
2336#ifdef COIN_HAS_CLP
2337    OsiClpSolverInterface *clpSolver
2338      = dynamic_cast<OsiClpSolverInterface *>(solver_);
2339    if (clpSolver) {
2340      int specialOptions = clpSolver->getModelPtr()->specialOptions();
2341      // 2097152 switches on rays in crunch
2342      if (!parentModel_)
2343        clpSolver->getModelPtr()->setSpecialOptions(specialOptions | 32 | 2097152);
2344      else
2345        clpSolver->getModelPtr()->setSpecialOptions(specialOptions & ~(32 | 2097152));
2346    }
2347  }
2348#endif
2349#endif
2350#ifdef COIN_HAS_NTY
2351  // maybe allow on fix and restart later
2352  if ((moreSpecialOptions2_ & (128 | 256)) != 0 && !parentModel_) {
2353    symmetryInfo_ = new CbcSymmetry();
2354    symmetryInfo_->setupSymmetry(*continuousSolver_);
2355    int numberGenerators = symmetryInfo_->statsOrbits(this, 0);
2356    if (!symmetryInfo_->numberUsefulOrbits() && (moreSpecialOptions2_ & (128 | 256)) != (128 | 256)) {
2357      delete symmetryInfo_;
2358      symmetryInfo_ = NULL;
2359      moreSpecialOptions2_ &= ~(128 | 256);
2360    }
2361    if ((moreSpecialOptions2_ & (128 | 256)) == (128 | 256)) {
2362      //moreSpecialOptions2_ &= ~256;
2363    }
2364  }
2365#endif
2366
2367  // add cutoff as constraint if wanted
2368  if (cutoffRowNumber_ == -2) {
2369    if (!parentModel_) {
2370      int numberColumns = solver_->getNumCols();
2371      double *obj = CoinCopyOfArray(solver_->getObjCoefficients(), numberColumns);
2372      int *indices = new int[numberColumns];
2373      int n = 0;
2374      for (int i = 0; i < numberColumns; i++) {
2375        if (obj[i]) {
2376          indices[n] = i;
2377          obj[n++] = obj[i];
2378        }
2379      }
2380      if (n) {
2381        double cutoff = getCutoff();
2382        // relax a little bit
2383        cutoff += 1.0e-4;
2384        double offset;
2385        solver_->getDblParam(OsiObjOffset, offset);
2386        cutoffRowNumber_ = solver_->getNumRows();
2387        solver_->addRow(n, indices, obj, -COIN_DBL_MAX, CoinMin(cutoff, 1.0e25) + offset);
2388      } else {
2389        // no objective!
2390        cutoffRowNumber_ = -1;
2391      }
2392      delete[] indices;
2393      delete[] obj;
2394    } else {
2395      // switch off
2396      cutoffRowNumber_ = -1;
2397    }
2398  }
2399  numberRowsAtContinuous_ = getNumRows();
2400  solver_->saveBaseModel();
2401  /*
2402      Check the objective to see if we can deduce a nontrivial increment. If
2403      it's better than the current value for CbcCutoffIncrement, it'll be
2404      installed.
2405    */
2406  if (solverCharacteristics_->reducedCostsAccurate())
2407    analyzeObjective();
2408  {
2409    // may be able to change cutoff now
2410    double cutoff = getCutoff();
2411    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
2412    if (cutoff > bestObjective_ - increment) {
2413      cutoff = bestObjective_ - increment;
2414      setCutoff(cutoff);
2415    }
2416  }
2417#ifdef COIN_HAS_CLP
2418  // Possible save of pivot method
2419  ClpDualRowPivot *savePivotMethod = NULL;
2420  {
2421    // pass tolerance and increment to solver
2422    OsiClpSolverInterface *clpSolver
2423      = dynamic_cast<OsiClpSolverInterface *>(solver_);
2424    if (clpSolver)
2425      clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2426#ifdef CLP_RESOLVE
2427    if ((moreSpecialOptions_ & 1048576) != 0 && !parentModel_ && clpSolver) {
2428      resolveClp(clpSolver, 0);
2429    }
2430#endif
2431  }
2432#endif
2433  /*
2434      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2435      the active subproblem. whichGenerator will be used to record the generator
2436      that produced a given cut.
2437    */
2438#define INITIAL_MAXIMUM_WHICH 1000
2439  maximumWhich_ = INITIAL_MAXIMUM_WHICH;
2440  delete[] whichGenerator_;
2441  whichGenerator_ = new int[maximumWhich_];
2442  memset(whichGenerator_, 0, maximumWhich_ * sizeof(int));
2443  maximumNumberCuts_ = 0;
2444  currentNumberCuts_ = 0;
2445  delete[] addedCuts_;
2446  addedCuts_ = NULL;
2447  OsiObject **saveObjects = NULL;
2448  maximumRows_ = numberRowsAtContinuous_;
2449  currentDepth_ = 0;
2450  workingBasis_.resize(maximumRows_, numberColumns);
2451  /*
2452      Set up an empty heap and associated data structures to hold the live set
2453      (problems which require further exploration).
2454    */
2455  CbcCompareDefault *compareActual
2456    = dynamic_cast<CbcCompareDefault *>(nodeCompare_);
2457  if (compareActual) {
2458    compareActual->setBestPossible(direction * solver_->getObjValue());
2459    compareActual->setCutoff(getCutoff());
2460#ifdef JJF_ZERO
2461    if (false && !numberThreads_ && !parentModel_) {
2462      printf("CbcTreeArray ? threads ? parentArray\n");
2463      // Setup new style tree
2464      delete tree_;
2465      tree_ = new CbcTreeArray();
2466    }
2467#endif
2468  }
2469  tree_->setComparison(*nodeCompare_);
2470  /*
2471      Used to record the path from a node to the root of the search tree, so that
2472      we can then traverse from the root to the node when restoring a subproblem.
2473    */
2474  maximumDepth_ = 10;
2475  delete[] walkback_;
2476  walkback_ = new CbcNodeInfo *[maximumDepth_];
2477  lastDepth_ = 0;
2478  delete[] lastNodeInfo_;
2479  lastNodeInfo_ = new CbcNodeInfo *[maximumDepth_];
2480  delete[] lastNumberCuts_;
2481  lastNumberCuts_ = new int[maximumDepth_];
2482  maximumCuts_ = 100;
2483  lastNumberCuts2_ = 0;
2484  delete[] lastCut_;
2485  lastCut_ = new const OsiRowCut *[maximumCuts_];
2486  /*
2487      Used to generate bound edits for CbcPartialNodeInfo.
2488    */
2489  double *lowerBefore = new double[numberColumns];
2490  double *upperBefore = new double[numberColumns];
2491  /*
2492    Set up to run heuristics and generate cuts at the root node. The heavy
2493    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
2494
2495    To start, tell cut generators they can be a bit more aggressive at the
2496    root node.
2497
2498    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2499        with cuts at root'. Is phase_ = 1 the correct indication when
2500        doHeurisiticsAtRoot is called to run heuristics outside of the main
2501        cut / heurisitc / reoptimise loop in solveWithCuts?
2502
2503      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2504      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2505      fixing) until no cuts are generated, the change in objective falls off,  or
2506      the limit on the number of rounds of cut generation is exceeded.
2507
2508      At the end of all this, any cuts will be recorded in cuts and also
2509      installed in the solver's constraint system. We'll have reoptimised, and
2510      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2511      adjusted accordingly).
2512
2513      Tell cut generators they can be a bit more aggressive at root node
2514
2515      TODO: Why don't we make a copy of the solution after solveWithCuts?
2516      TODO: If numberUnsatisfied == 0, don't we have a solution?
2517    */
2518  phase_ = 1;
2519  int iCutGenerator;
2520  for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2521    // If parallel switch off global cuts
2522    if (numberThreads_) {
2523      generator_[iCutGenerator]->setGlobalCuts(false);
2524      generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2525    }
2526    CglCutGenerator *generator = generator_[iCutGenerator]->generator();
2527    generator->setAggressiveness(generator->getAggressiveness() + 100);
2528    if (!generator->canDoGlobalCuts())
2529      generator->setGlobalCuts(false);
2530  }
2531  OsiCuts cuts;
2532  int anyAction = -1;
2533  numberOldActiveCuts_ = 0;
2534  numberNewCuts_ = 0;
2535  // Array to mark solution
2536  delete[] usedInSolution_;
2537  usedInSolution_ = new int[numberColumns];
2538  CoinZeroN(usedInSolution_, numberColumns);
2539  /*
2540      For printing totals and for CbcNode (numberNodes_)
2541    */
2542  numberIterations_ = 0;
2543  numberNodes_ = 0;
2544  numberNodes2_ = 0;
2545  maximumStatistics_ = 0;
2546  maximumDepthActual_ = 0;
2547  numberDJFixed_ = 0.0;
2548  if (!parentModel_) {
2549    if ((specialOptions_ & 262144) != 0) {
2550      // create empty stored cuts
2551      //storedRowCuts_ = new CglStored(solver_->getNumCols());
2552    } else if ((specialOptions_ & 524288) != 0 && storedRowCuts_) {
2553      // tighten and set best solution
2554      // A) tight bounds on integer variables
2555      /*
2556                storedRowCuts_ are coming in from outside, probably for nonlinear.
2557              John was unsure about origin.
2558            */
2559      const double *lower = solver_->getColLower();
2560      const double *upper = solver_->getColUpper();
2561      const double *tightLower = storedRowCuts_->tightLower();
2562      const double *tightUpper = storedRowCuts_->tightUpper();
2563      int nTightened = 0;
2564      for (int i = 0; i < numberIntegers_; i++) {
2565        int iColumn = integerVariable_[i];
2566        if (tightLower[iColumn] > lower[iColumn]) {
2567          nTightened++;
2568          solver_->setColLower(iColumn, tightLower[iColumn]);
2569        }
2570        if (tightUpper[iColumn] < upper[iColumn]) {
2571          nTightened++;
2572          solver_->setColUpper(iColumn, tightUpper[iColumn]);
2573        }
2574      }
2575      if (nTightened)
2576        COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2577      if (storedRowCuts_->bestObjective() < bestObjective_) {
2578        // B) best solution
2579        double objValue = storedRowCuts_->bestObjective();
2580        setBestSolution(CBC_SOLUTION, objValue,
2581          storedRowCuts_->bestSolution());
2582        // Do heuristics
2583        // Allow RINS
2584        for (int i = 0; i < numberHeuristics_; i++) {
2585          CbcHeuristicRINS *rins
2586            = dynamic_cast<CbcHeuristicRINS *>(heuristic_[i]);
2587          if (rins) {
2588            rins->setLastNode(-100);
2589          }
2590        }
2591      }
2592    }
2593  }
2594#ifdef SWITCH_VARIABLES
2595  // see if any switching variables
2596  if (numberIntegers_ < solver_->getNumCols())
2597    findSwitching();
2598#endif
2599  /*
2600      Run heuristics at the root. This is the only opportunity to run FPump; it
2601      will be removed from the heuristics list by doHeuristicsAtRoot.
2602    */
2603  // See if multiple runs wanted
2604  CbcModel **rootModels = NULL;
2605  if (!parentModel_ && multipleRootTries_ % 100) {
2606    double rootTimeCpu = CoinCpuTime();
2607    double startTimeRoot = CoinGetTimeOfDay();
2608    int numberRootThreads = 1;
2609    /* undocumented fine tuning
2610         aabbcc where cc is number of tries
2611         bb if nonzero is number of threads
2612         aa if nonzero just do heuristics
2613      */
2614    int numberModels = multipleRootTries_ % 100;
2615#ifdef CBC_THREAD
2616    numberRootThreads = (multipleRootTries_ / 100) % 100;
2617    if (!numberRootThreads) {
2618      if (numberThreads_ < 2)
2619        numberRootThreads = numberModels;
2620      else
2621        numberRootThreads = CoinMin(numberThreads_, numberModels);
2622    }
2623#endif
2624    int otherOptions = (multipleRootTries_ / 10000) % 100;
2625    rootModels = new CbcModel *[numberModels];
2626    unsigned int newSeed = randomSeed_;
2627    if (newSeed == 0) {
2628      double time = fabs(CoinGetTimeOfDay());
2629      while (time >= COIN_INT_MAX)
2630        time *= 0.5;
2631      newSeed = static_cast<unsigned int>(time);
2632    } else if (newSeed < 0) {
2633      newSeed = 123456789;
2634#ifdef COIN_HAS_CLP
2635      OsiClpSolverInterface *clpSolver
2636        = dynamic_cast<OsiClpSolverInterface *>(solver_);
2637      if (clpSolver) {
2638        newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2639      }
2640#endif
2641    }
2642    CoinWarmStartBasis *basis = dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart());
2643    for (int i = 0; i < numberModels; i++) {
2644      rootModels[i] = new CbcModel(*this);
2645      rootModels[i]->setNumberThreads(0);
2646      rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2647      rootModels[i]->setRandomSeed(newSeed + 10000000 * i);
2648      rootModels[i]->randomNumberGenerator()->setSeed(newSeed + 50000000 * i);
2649      rootModels[i]->setMultipleRootTries(0);
2650#ifdef COIN_HAS_NTY
2651      rootModels[i]->zapSymmetry();
2652      rootModels[i]->moreSpecialOptions2_ &= ~(128 | 256); // off nauty
2653#endif
2654      // use seed
2655      rootModels[i]->setSpecialOptions(specialOptions_ | (4194304 | 8388608));
2656      rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & (~(134217728 | 4194304)));
2657      rootModels[i]->setMoreSpecialOptions2(moreSpecialOptions2_ & (~(128 | 256)));
2658      rootModels[i]->solver_->setWarmStart(basis);
2659#ifdef COIN_HAS_CLP
2660      OsiClpSolverInterface *clpSolver
2661        = dynamic_cast<OsiClpSolverInterface *>(rootModels[i]->solver_);
2662#define NEW_RANDOM_BASIS
2663#ifdef NEW_RANDOM_BASIS
2664      if (i == 0)
2665        continue;
2666#endif
2667      if (clpSolver) {
2668        ClpSimplex *simplex = clpSolver->getModelPtr();
2669        if (defaultHandler_)
2670          simplex->setDefaultMessageHandler();
2671        simplex->setRandomSeed(newSeed + 20000000 * i);
2672        simplex->allSlackBasis();
2673        int logLevel = simplex->logLevel();
2674        if (logLevel == 1)
2675          simplex->setLogLevel(0);
2676        if (i != 0) {
2677#ifdef NEW_RANDOM_BASIS
2678          int numberRows = simplex->numberRows();
2679          int throwOut = 20; //2+numberRows/100;
2680          for (int iThrow = 0; iThrow < throwOut; iThrow++) {
2681            double random = simplex->randomNumberGenerator()->randomDouble();
2682            int iStart = static_cast<int>(random * numberRows);
2683            for (int j = iStart; j < numberRows; j++) {
2684              if (simplex->getRowStatus(j) != ClpSimplex::basic) {
2685                simplex->setRowStatus(j, ClpSimplex::basic);
2686                break;
2687              }
2688            }
2689          }
2690          clpSolver->setWarmStart(NULL);
2691#else
2692            double random = simplex->randomNumberGenerator()->randomDouble();
2693            int bias = static_cast<int>(random * (numberIterations / 4));
2694            simplex->setMaximumIterations(numberIterations / 2 + bias);
2695            simplex->primal();
2696            simplex->setMaximumIterations(COIN_INT_MAX);
2697            simplex->dual();
2698#endif
2699        } else {
2700#ifndef NEW_RANDOM_BASIS
2701          simplex->primal();
2702#endif
2703#endif
2704        }
2705#ifdef NEW_RANDOM_BASIS
2706        simplex->setLogLevel(logLevel);
2707        clpSolver->setWarmStart(NULL);
2708#endif
2709      }
2710      for (int j = 0; j < numberHeuristics_; j++)
2711        rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed() + 100000000 * i);
2712      for (int j = 0; j < numberCutGenerators_; j++)
2713        rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
2714    }
2715    delete basis;
2716#ifdef CBC_THREAD
2717    if (numberRootThreads == 1) {
2718#endif
2719      for (int iModel = 0; iModel < numberModels; iModel++) {
2720        doRootCbcThread(rootModels[iModel]);
2721        // see if solved at root node
2722        if (rootModels[iModel]->getMaximumNodes()) {
2723          feasible = false;
2724          break;
2725        }
2726      }
2727#ifdef CBC_THREAD
2728    } else {
2729      Coin_pthread_t *threadId = new Coin_pthread_t[numberRootThreads];
2730      for (int kModel = 0; kModel < numberModels; kModel += numberRootThreads) {
2731        bool finished = false;
2732        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
2733          pthread_create(&(threadId[iModel - kModel].thr), NULL,
2734            doRootCbcThread,
2735            rootModels[iModel]);
2736        }
2737        // wait
2738        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
2739          pthread_join(threadId[iModel - kModel].thr, NULL);
2740        }
2741        // see if solved at root node
2742        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
2743          if (rootModels[iModel]->getMaximumNodes())
2744            finished = true;
2745        }
2746        if (finished) {
2747          feasible = false;
2748          break;
2749        }
2750      }
2751      delete[] threadId;
2752    }
2753#endif
2754    // sort solutions
2755    int *which = new int[numberModels];
2756    double *value = new double[numberModels];
2757    int numberSolutions = 0;
2758    for (int iModel = 0; iModel < numberModels; iModel++) {
2759      if (rootModels[iModel]->bestSolution()) {
2760        which[numberSolutions] = iModel;
2761        value[numberSolutions++] = -rootModels[iModel]->getMinimizationObjValue();
2762      }
2763    }
2764    char general[100];
2765    rootTimeCpu = CoinCpuTime() - rootTimeCpu;
2766    if (numberRootThreads == 1)
2767      sprintf(general, "Multiple root solvers took a total of %.2f seconds\n",
2768        rootTimeCpu);
2769    else
2770      sprintf(general, "Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2771        rootTimeCpu, CoinGetTimeOfDay() - startTimeRoot);
2772    messageHandler()->message(CBC_GENERAL,
2773      messages())
2774      << general << CoinMessageEol;
2775    CoinSort_2(value, value + numberSolutions, which);
2776    // to get name
2777    CbcHeuristicRINS dummyHeuristic;
2778    dummyHeuristic.setHeuristicName("Multiple root solvers");
2779    lastHeuristic_ = &dummyHeuristic;
2780    for (int i = 0; i < numberSolutions; i++) {
2781      double objValue = -value[i];
2782      if (objValue < getCutoff()) {
2783        int iModel = which[i];
2784        setBestSolution(CBC_ROUNDING, objValue,
2785          rootModels[iModel]->bestSolution());
2786      }
2787    }
2788    lastHeuristic_ = NULL;
2789    delete[] which;
2790    delete[] value;
2791  }
2792  // Do heuristics
2793  if (numberObjects_ && !rootModels)
2794    doHeuristicsAtRoot();
2795  if (solverCharacteristics_->solutionAddsCuts()) {
2796    // With some heuristics solver needs a resolve here
2797    solver_->resolve();
2798    if (!isProvenOptimal()) {
2799      solver_->initialSolve();
2800    }
2801  }
2802  /*
2803      Grepping through the code, it would appear that this is a command line
2804      debugging hook.  There's no obvious place in the code where this is set to
2805      a negative value.
2806
2807      User hook, says John.
2808    */
2809  if (intParam_[CbcMaxNumNode] < 0
2810    || numberSolutions_ >= getMaximumSolutions())
2811    eventHappened_ = true; // stop as fast as possible
2812  stoppedOnGap_ = false;
2813  // See if can stop on gap
2814  bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2815  if (canStopOnGap()) {
2816    if (bestPossibleObjective_ < getCutoff())
2817      stoppedOnGap_ = true;
2818    feasible = false;
2819    //eventHappened_=true; // stop as fast as possible
2820  }
2821  /*
2822      Set up for statistics collection, if requested. Standard values are
2823      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2824      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2825      command line debugging hook.
2826    */
2827  statistics_ = NULL;
2828  // Do on switch
2829  if (doStatistics > 0 && doStatistics <= 100) {
2830    maximumStatistics_ = 10000;
2831    statistics_ = new CbcStatistics *[maximumStatistics_];
2832    memset(statistics_, 0, maximumStatistics_ * sizeof(CbcStatistics *));
2833  }
2834  // See if we can add integers
2835  if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_ & 65536) != 0 && !parentModel_ && false) {
2836    int numberIntegers1 = 0;
2837    int numberColumns = solver_->getNumCols();
2838    for (int i = 0; i < numberColumns; i++) {
2839      if (solver_->isInteger(i))
2840        numberIntegers1++;
2841    }
2842    AddIntegers();
2843    // make sure in sync
2844    int numberIntegers2 = 0;
2845    for (int i = 0; i < numberColumns; i++) {
2846      if (solver_->isInteger(i))
2847        numberIntegers2++;
2848    }
2849    if (numberIntegers1 < numberIntegers2) {
2850      findIntegers(true, 2);
2851      convertToDynamic();
2852    }
2853  }
2854
2855  /*
2856      Do an initial round of cut generation for the root node. Depending on the
2857      type of underlying solver, we may want to do this even if the initial query
2858      to the objects indicates they're satisfied.
2859
2860      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2861      loop (including reduced cost fixing) until no cuts are generated, the
2862      change in objective falls off,  or the limit on the number of rounds of cut
2863      generation is exceeded.
2864
2865      At the end of all this, any cuts will be recorded in cuts and also
2866      installed in the solver's constraint system. We'll have reoptimised, and
2867      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2868      adjusted accordingly).
2869    */
2870  int iObject;
2871  int numberUnsatisfied = 0;
2872  delete[] currentSolution_;
2873  currentSolution_ = new double[numberColumns];
2874  testSolution_ = currentSolution_;
2875  memcpy(currentSolution_, solver_->getColSolution(),
2876    numberColumns * sizeof(double));
2877  // point to useful information
2878  OsiBranchingInformation usefulInfo = usefulInformation();
2879
2880  for (iObject = 0; iObject < numberObjects_; iObject++) {
2881    double infeasibility = object_[iObject]->checkInfeasibility(&usefulInfo);
2882    if (infeasibility)
2883      numberUnsatisfied++;
2884  }
2885  // replace solverType
2886  double *tightBounds = NULL;
2887  if (solverCharacteristics_->tryCuts()) {
2888
2889    if (numberUnsatisfied) {
2890      // User event
2891      if (!eventHappened_ && feasible) {
2892        if (rootModels) {
2893          // for fixings
2894          int numberColumns = solver_->getNumCols();
2895          tightBounds = new double[2 * numberColumns];
2896          {
2897            const double *lower = solver_->getColLower();
2898            const double *upper = solver_->getColUpper();
2899            for (int i = 0; i < numberColumns; i++) {
2900              tightBounds[2 * i + 0] = lower[i];
2901              tightBounds[2 * i + 1] = upper[i];
2902            }
2903          }
2904          int numberModels = multipleRootTries_ % 100;
2905          const OsiSolverInterface **solvers = new const OsiSolverInterface *[numberModels];
2906          int numberRows = continuousSolver_->getNumRows();
2907          int maxCuts = 0;
2908          for (int i = 0; i < numberModels; i++) {
2909            solvers[i] = rootModels[i]->solver();
2910            const double *lower = solvers[i]->getColLower();
2911            const double *upper = solvers[i]->getColUpper();
2912            for (int j = 0; j < numberColumns; j++) {
2913              tightBounds[2 * j + 0] = CoinMax(lower[j], tightBounds[2 * j + 0]);
2914              tightBounds[2 * j + 1] = CoinMin(upper[j], tightBounds[2 * j + 1]);
2915            }
2916            int numberRows2 = solvers[i]->getNumRows();
2917            assert(numberRows2 >= numberRows);
2918            maxCuts += numberRows2 - numberRows;
2919            // accumulate statistics
2920            for (int j = 0; j < numberCutGenerators_; j++) {
2921              generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2922            }
2923          }
2924          for (int j = 0; j < numberCutGenerators_; j++) {
2925            generator_[j]->scaleBackStatistics(numberModels);
2926          }
2927          //CbcRowCuts rowCut(maxCuts);
2928          const OsiRowCutDebugger *debugger = NULL;
2929          if ((specialOptions_ & 1) != 0)
2930            debugger = solver_->getRowCutDebugger();
2931          for (int iModel = 0; iModel < numberModels; iModel++) {
2932            int numberRows2 = solvers[iModel]->getNumRows();
2933            const CoinPackedMatrix *rowCopy = solvers[iModel]->getMatrixByRow();
2934            const int *rowLength = rowCopy->getVectorLengths();
2935            const double *elements = rowCopy->getElements();
2936            const int *column = rowCopy->getIndices();
2937            const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
2938            const double *rowLower = solvers[iModel]->getRowLower();
2939            const double *rowUpper = solvers[iModel]->getRowUpper();
2940            for (int iRow = numberRows; iRow < numberRows2; iRow++) {
2941              OsiRowCut rc;
2942              rc.setLb(rowLower[iRow]);
2943              rc.setUb(rowUpper[iRow]);
2944              CoinBigIndex start = rowStart[iRow];
2945              rc.setRow(rowLength[iRow], column + start, elements + start, false);
2946              if (debugger)
2947                CoinAssert(!debugger->invalidCut(rc));
2948              globalCuts_.addCutIfNotDuplicate(rc);
2949            }
2950            //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
2951            //numberCuts += cutsAdded;
2952            //printf("Model %d gave %d cuts (out of %d possible)\n",
2953            //     iModel,cutsAdded,numberRows2-numberRows);
2954          }
2955          // normally replace global cuts
2956          //if (!globalCuts_.())
2957          //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2958          //rowCut.addCuts(globalCuts_);
2959          int nTightened = 0;
2960          assert(feasible);
2961          {
2962            double tolerance = 1.0e-5;
2963            const double *lower = solver_->getColLower();
2964            const double *upper = solver_->getColUpper();
2965            for (int i = 0; i < numberColumns; i++) {
2966              if (tightBounds[2 * i + 0] > tightBounds[2 * i + 1] + 1.0e-9) {
2967                feasible = false;
2968                char general[200];
2969                sprintf(general, "Solvers give infeasible bounds on %d %g,%g was %g,%g - search finished\n",
2970                  i, tightBounds[2 * i + 0], tightBounds[2 * i + 1], lower[i], upper[i]);
2971                messageHandler()->message(CBC_GENERAL, messages())
2972                  << general << CoinMessageEol;
2973                break;
2974              }
2975              double oldLower = lower[i];
2976              double oldUpper = upper[i];
2977              if (tightBounds[2 * i + 0] > oldLower + tolerance) {
2978                nTightened++;
2979                solver_->setColLower(i, tightBounds[2 * i + 0]);
2980              }
2981              if (tightBounds[2 * i + 1] < oldUpper - tolerance) {
2982                nTightened++;
2983                solver_->setColUpper(i, tightBounds[2 * i + 1]);
2984              }
2985            }
2986          }
2987          delete[] tightBounds;
2988          tightBounds = NULL;
2989          char printBuffer[200];
2990          sprintf(printBuffer, "%d solvers added %d different cuts out of pool of %d",
2991            numberModels, globalCuts_.sizeRowCuts(), maxCuts);
2992          messageHandler()->message(CBC_GENERAL, messages())
2993            << printBuffer << CoinMessageEol;
2994          if (nTightened) {
2995            sprintf(printBuffer, "%d bounds were tightened",
2996              nTightened);
2997            messageHandler()->message(CBC_GENERAL, messages())
2998              << printBuffer << CoinMessageEol;
2999          }
3000          delete[] solvers;
3001        }
3002        if (!parentModel_ && (moreSpecialOptions_ & 67108864) != 0) {
3003          // load cuts from file
3004          FILE *fp = fopen("global.cuts", "rb");
3005          if (fp) {
3006            size_t nRead;
3007            int numberColumns = solver_->getNumCols();
3008            int numCols;
3009            nRead = fread(&numCols, sizeof(int), 1, fp);
3010            if (nRead != 1)
3011              throw("Error in fread");
3012            if (numberColumns != numCols) {
3013              printf("Mismatch on columns %d %d\n", numberColumns, numCols);
3014              fclose(fp);
3015            } else {
3016              // If rootModel just do some
3017              double threshold = -1.0;
3018              if (!multipleRootTries_)
3019                threshold = 0.5;
3020              int initialCuts = 0;
3021              int initialGlobal = globalCuts_.sizeRowCuts();
3022              double *elements = new double[numberColumns + 2];
3023              int *indices = new int[numberColumns];
3024              int numberEntries = 1;
3025              while (numberEntries > 0) {
3026                nRead = fread(&numberEntries, sizeof(int), 1, fp);
3027                if (nRead != 1)
3028                  throw("Error in fread");
3029                double randomNumber = randomNumberGenerator_.randomDouble();
3030                if (numberEntries > 0) {
3031                  initialCuts++;
3032                  nRead = fread(elements, sizeof(double), numberEntries + 2, fp);
3033                  if (nRead != static_cast<size_t>(numberEntries + 2))
3034                    throw("Error in fread");
3035                  nRead = fread(indices, sizeof(int), numberEntries, fp);
3036                  if (nRead != static_cast<size_t>(numberEntries))
3037                    throw("Error in fread");
3038                  if (randomNumber > threshold) {
3039                    OsiRowCut rc;
3040                    rc.setLb(elements[numberEntries]);
3041                    rc.setUb(elements[numberEntries + 1]);
3042                    rc.setRow(numberEntries, indices, elements,
3043                      false);
3044                    rc.setGloballyValidAsInteger(2);
3045                    globalCuts_.addCutIfNotDuplicate(rc);
3046                  }
3047                }
3048              }
3049              fclose(fp);
3050              // fixes
3051              int nTightened = 0;
3052              fp = fopen("global.fix", "rb");
3053              if (fp) {
3054                nRead = fread(indices, sizeof(int), 2, fp);
3055                if (nRead != 2)
3056                  throw("Error in fread");
3057                if (numberColumns != indices[0]) {
3058                  printf("Mismatch on columns %d %d\n", numberColumns,
3059                    indices[0]);
3060                } else {
3061                  indices[0] = 1;
3062                  while (indices[0] >= 0) {
3063                    nRead = fread(indices, sizeof(int), 2, fp);
3064                    if (nRead != 2)
3065                      throw("Error in fread");
3066                    int iColumn = indices[0];
3067                    if (iColumn >= 0) {
3068                      nTightened++;
3069                      nRead = fread(elements, sizeof(double), 4, fp);
3070                      if (nRead != 4)
3071                        throw("Error in fread");
3072                      solver_->setColLower(iColumn, elements[0]);
3073                      solver_->setColUpper(iColumn, elements[1]);
3074                    }
3075                  }
3076                }
3077              }
3078              if (fp)
3079                fclose(fp);
3080              char printBuffer[200];
3081              sprintf(printBuffer, "%d cuts read in of which %d were unique, %d bounds tightened",
3082                initialCuts,
3083                globalCuts_.sizeRowCuts() - initialGlobal, nTightened);
3084              messageHandler()->message(CBC_GENERAL, messages())
3085                << printBuffer << CoinMessageEol;
3086              delete[] elements;
3087              delete[] indices;
3088            }
3089          }
3090        }
3091        if (feasible)
3092          feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3093            NULL);
3094        if (multipleRootTries_ && (moreSpecialOptions_ & 134217728) != 0) {
3095          FILE *fp = NULL;
3096          size_t nRead;
3097          int numberColumns = solver_->getNumCols();
3098          int initialCuts = 0;
3099          if ((moreSpecialOptions_ & 134217728) != 0) {
3100            // append so go down to end
3101            fp = fopen("global.cuts", "r+b");
3102            if (fp) {
3103              int numCols;
3104              nRead = fread(&numCols, sizeof(int), 1, fp);
3105              if (nRead != 1)
3106                throw("Error in fread");
3107              if (numberColumns != numCols) {
3108                printf("Mismatch on columns %d %d\n", numberColumns, numCols);
3109                fclose(fp);
3110                fp = NULL;
3111              }
3112            }
3113          }
3114          double *elements = new double[numberColumns + 2];
3115          int *indices = new int[numberColumns];
3116          if (fp) {
3117            int numberEntries = 1;
3118            while (numberEntries > 0) {
3119              fpos_t position;
3120              fgetpos(fp, &position);
3121              nRead = fread(&numberEntries, sizeof(int), 1, fp);
3122              if (nRead != 1)
3123                throw("Error in fread");
3124              if (numberEntries > 0) {
3125                initialCuts++;
3126                nRead = fread(elements, sizeof(double), numberEntries + 2, fp);
3127                if (nRead != static_cast<size_t>(numberEntries + 2))
3128                  throw("Error in fread");
3129                nRead = fread(indices, sizeof(int), numberEntries, fp);
3130                if (nRead != static_cast<size_t>(numberEntries))
3131                  throw("Error in fread");
3132              } else {
3133                // end
3134                fsetpos(fp, &position);
3135              }
3136            }
3137          } else {
3138            fp = fopen("global.cuts", "wb");
3139            size_t nWrite;
3140            nWrite = fwrite(&numberColumns, sizeof(int), 1, fp);
3141            if (nWrite != 1)
3142              throw("Error in fwrite");
3143          }
3144          size_t nWrite;
3145          // now append binding cuts
3146          int numberC = continuousSolver_->getNumRows();
3147          int numberRows = solver_->getNumRows();
3148          printf("Saving %d cuts (up from %d)\n",
3149            initialCuts + numberRows - numberC, initialCuts);
3150          const double *rowLower = solver_->getRowLower();
3151          const double *rowUpper = solver_->getRowUpper();
3152          // Row copy
3153          CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3154          const double *elementByRow = matrixByRow.getElements();
3155          const int *column = matrixByRow.getIndices();
3156          const CoinBigIndex *rowStart = matrixByRow.getVectorStarts();
3157          const int *rowLength = matrixByRow.getVectorLengths();
3158          for (int iRow = numberC; iRow < numberRows; iRow++) {
3159            int n = rowLength[iRow];
3160            assert(n);
3161            CoinBigIndex start = rowStart[iRow];
3162            memcpy(elements, elementByRow + start, n * sizeof(double));
3163            memcpy(indices, column + start, n * sizeof(int));
3164            elements[n] = rowLower[iRow];
3165            elements[n + 1] = rowUpper[iRow];
3166            nWrite = fwrite(&n, sizeof(int), 1, fp);
3167            if (nWrite != 1)
3168              throw("Error in fwrite");
3169            nWrite = fwrite(elements, sizeof(double), n + 2, fp);
3170            if (nWrite != static_cast<size_t>(n + 2))
3171              throw("Error in fwrite");
3172            nWrite = fwrite(indices, sizeof(int), n, fp);
3173            if (nWrite != static_cast<size_t>(n))
3174              throw("Error in fwrite");
3175          }
3176          // eof marker
3177          int eofMarker = -1;
3178          nWrite = fwrite(&eofMarker, sizeof(int), 1, fp);
3179          if (nWrite != 1)
3180            throw("Error in fwrite");
3181          fclose(fp);
3182          // do tighter bounds (? later extra to original columns)
3183          int nTightened = 0;
3184          const double *lower = solver_->getColLower();
3185          const double *upper = solver_->getColUpper();
3186          const double *originalLower = continuousSolver_->getColLower();
3187          const double *originalUpper = continuousSolver_->getColUpper();
3188          double tolerance = 1.0e-5;
3189          for (int i = 0; i < numberColumns; i++) {
3190            if (lower[i] > originalLower[i] + tolerance) {
3191              nTightened++;
3192            }
3193            if (upper[i] < originalUpper[i] - tolerance) {
3194              nTightened++;
3195            }
3196          }
3197          if (nTightened) {
3198            fp = fopen("global.fix", "wb");
3199            size_t nWrite;
3200            indices[0] = numberColumns;
3201            if (originalColumns_)
3202              indices[1] = COIN_INT_MAX;
3203            else
3204              indices[1] = -1;
3205            nWrite = fwrite(indices, sizeof(int), 2, fp);
3206            if (nWrite != 2)
3207              throw("Error in fwrite");
3208            for (int i = 0; i < numberColumns; i++) {
3209              int nTightened = 0;
3210              if (lower[i] > originalLower[i] + tolerance) {
3211                nTightened++;
3212              }
3213              if (upper[i] < originalUpper[i] - tolerance) {
3214                nTightened++;
3215              }
3216              if (nTightened) {
3217                indices[0] = i;
3218                if (originalColumns_)
3219                  indices[1] = originalColumns_[i];
3220                elements[0] = lower[i];
3221                elements[1] = upper[i];
3222                elements[2] = originalLower[i];
3223                elements[3] = originalUpper[i];
3224                nWrite = fwrite(indices, sizeof(int), 2, fp);
3225                if (nWrite != 2)
3226                  throw("Error in fwrite");
3227                nWrite = fwrite(elements, sizeof(double), 4, fp);
3228                if (nWrite != 4)
3229                  throw("Error in fwrite");
3230              }
3231            }
3232            // eof marker
3233            indices[0] = -1;
3234            nWrite = fwrite(indices, sizeof(int), 2, fp);
3235            if (nWrite != 2)
3236              throw("Error in fwrite");
3237            fclose(fp);
3238          }
3239          delete[] elements;
3240          delete[] indices;
3241        }
3242        if ((specialOptions_ & 524288) != 0 && !parentModel_
3243          && storedRowCuts_) {
3244          if (feasible) {
3245            /* pick up stuff and try again
3246                        add cuts, maybe keep around
3247                        do best solution and if so new heuristics
3248                        obviously tighten bounds
3249                        */
3250            // A and B probably done on entry
3251            // A) tight bounds on integer variables
3252            const double *lower = solver_->getColLower();
3253            const double *upper = solver_->getColUpper();
3254            const double *tightLower = storedRowCuts_->tightLower();
3255            const double *tightUpper = storedRowCuts_->tightUpper();
3256            int nTightened = 0;
3257            for (int i = 0; i < numberIntegers_; i++) {
3258              int iColumn = integerVariable_[i];
3259              if (tightLower[iColumn] > lower[iColumn]) {
3260                nTightened++;
3261                solver_->setColLower(iColumn, tightLower[iColumn]);
3262              }
3263              if (tightUpper[iColumn] < upper[iColumn]) {
3264                nTightened++;
3265                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3266              }
3267            }
3268            if (nTightened)
3269              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3270            if (storedRowCuts_->bestObjective() < bestObjective_) {
3271              // B) best solution
3272              double objValue = storedRowCuts_->bestObjective();
3273              setBestSolution(CBC_SOLUTION, objValue,
3274                storedRowCuts_->bestSolution());
3275              // Do heuristics
3276              // Allow RINS
3277              for (int i = 0; i < numberHeuristics_; i++) {
3278                CbcHeuristicRINS *rins
3279                  = dynamic_cast<CbcHeuristicRINS *>(heuristic_[i]);
3280                if (rins) {
3281                  rins->setLastNode(-100);
3282                }
3283              }
3284              doHeuristicsAtRoot();
3285            }
3286#ifdef JJF_ZERO
3287            int nCuts = storedRowCuts_->sizeRowCuts();
3288            // add to global list
3289            for (int i = 0; i < nCuts; i++) {
3290              OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3291              newCut.setGloballyValidAsInteger(2);
3292              newCut.mutableRow().setTestForDuplicateIndex(false);
3293              globalCuts_.insert(newCut);
3294            }
3295#else
3296          addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3297            true, false, false, -200);
3298#endif
3299            // Set cuts as active
3300            delete[] addedCuts_;
3301            maximumNumberCuts_ = cuts.sizeRowCuts();
3302            if (maximumNumberCuts_) {
3303              addedCuts_ = new CbcCountRowCut *[maximumNumberCuts_];
3304            } else {
3305              addedCuts_ = NULL;
3306            }
3307            for (int i = 0; i < maximumNumberCuts_; i++)
3308              addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3309                NULL, -1, -1, 2);
3310            COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3311            cuts = OsiCuts();
3312            currentNumberCuts_ = maximumNumberCuts_;
3313            feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3314              NULL);
3315            for (int i = 0; i < maximumNumberCuts_; i++)
3316              delete addedCuts_[i];
3317          }
3318          delete storedRowCuts_;
3319          storedRowCuts_ = NULL;
3320        }
3321      } else {
3322        feasible = false;
3323      }
3324    } else if (solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3325      // may generate cuts and turn the solution
3326      //to an infeasible one
3327      feasible = solveWithCuts(cuts, 2,
3328        NULL);
3329    }
3330  }
3331  if (rootModels) {
3332    int numberModels = multipleRootTries_ % 100;
3333    for (int i = 0; i < numberModels; i++)
3334      delete rootModels[i];
3335    delete[] rootModels;
3336  }
3337  // check extra info on feasibility
3338  if (!solverCharacteristics_->mipFeasible())
3339    feasible = false;
3340  // If max nodes==0 - don't do strong branching
3341  if (!getMaximumNodes()) {
3342    if (feasible)
3343      feasible = false;
3344    else
3345      setMaximumNodes(1); //allow to stop on success
3346  }
3347  topOfTree_ = NULL;
3348#ifdef CLP_RESOLVE
3349  if ((moreSpecialOptions_ & 2097152) != 0 && !parentModel_ && feasible) {
3350    OsiClpSolverInterface *clpSolver
3351      = dynamic_cast<OsiClpSolverInterface *>(solver_);
3352    if (clpSolver)
3353      resolveClp(clpSolver, 0);
3354  }
3355#endif
3356  // make cut generators less aggressive
3357  for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3358    CglCutGenerator *generator = generator_[iCutGenerator]->generator();
3359    generator->setAggressiveness(generator->getAggressiveness() - 100);
3360  }
3361  currentNumberCuts_ = numberNewCuts_;
3362  if (solverCharacteristics_->solutionAddsCuts()) {
3363    // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3364    solver_->resolve();
3365    if (!isProvenOptimal()) {
3366      solver_->initialSolve();
3367    }
3368  }
3369  // See if can stop on gap
3370  bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3371  if (canStopOnGap()) {
3372    if (bestPossibleObjective_ < getCutoff())
3373      stoppedOnGap_ = true;
3374    feasible = false;
3375  }
3376  // User event
3377  if (eventHappened_)
3378    feasible = false;
3379#if defined(COIN_HAS_CLP) && defined(COIN_HAS_CPX)
3380  /*
3381      This is the notion of using Cbc stuff to get going, then calling cplex to
3382      finish off.
3383    */
3384  if (feasible && (specialOptions_ & 16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3385    // Use Cplex to do search!
3386    double time1 = CoinCpuTime();
3387    OsiClpSolverInterface *clpSolver
3388      = dynamic_cast<OsiClpSolverInterface *>(solver_);
3389    OsiCpxSolverInterface cpxSolver;
3390    double direction = clpSolver->getObjSense();
3391    cpxSolver.setObjSense(direction);
3392    // load up cplex
3393    const CoinPackedMatrix *matrix = continuousSolver_->getMatrixByCol();
3394    const double *rowLower = continuousSolver_->getRowLower();
3395    const double *rowUpper = continuousSolver_->getRowUpper();
3396    const double *columnLower = continuousSolver_->getColLower();
3397    const double *columnUpper = continuousSolver_->getColUpper();
3398    const double *objective = continuousSolver_->getObjCoefficients();
3399    cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3400      objective, rowLower, rowUpper);
3401    double *setSol = new double[numberIntegers_];
3402    int *setVar = new int[numberIntegers_];
3403    // cplex doesn't know about objective offset
3404    double offset = clpSolver->getModelPtr()->objectiveOffset();
3405    for (int i = 0; i < numberIntegers_; i++) {
3406      int iColumn = integerVariable_[i];
3407      cpxSolver.setInteger(iColumn);
3408      if (bestSolution_) {
3409        setSol[i] = bestSolution_[iColumn];
3410        setVar[i] = iColumn;
3411      }
3412    }
3413    CPXENVptr env = cpxSolver.getEnvironmentPtr();
3414    CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3415    cpxSolver.switchToMIP();
3416    if (bestSolution_) {
3417#if 0
3418            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3419#else
3420        int zero = 0;
3421        CPXaddmipstarts(env, lpPtr, 1, numberIntegers_, &zero, setVar, setSol, NULL, NULL);
3422#endif
3423    }
3424    if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3425      // add cuts
3426      const CoinPackedMatrix *matrix = clpSolver->getMatrixByRow();
3427      const double *rhs = clpSolver->getRightHandSide();
3428      const char *rowSense = clpSolver->getRowSense();
3429      const double *elementByRow = matrix->getElements();
3430      const int *column = matrix->getIndices();
3431      const CoinBigIndex *rowStart = matrix->getVectorStarts();
3432      const int *rowLength = matrix->getVectorLengths();
3433      int nStart = continuousSolver_->getNumRows();
3434      int nRows = clpSolver->getNumRows();
3435      int size = rowStart[nRows - 1] + rowLength[nRows - 1] - rowStart[nStart];
3436      int nAdd = 0;
3437      double *rmatval = new double[size];
3438      int *rmatind = new int[size];
3439      int *rmatbeg = new int[nRows - nStart + 1];
3440      size = 0;
3441      rmatbeg[0] = 0;
3442      for (int i = nStart; i < nRows; i++) {
3443        for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3444          rmatind[size] = column[k];
3445          rmatval[size++] = elementByRow[k];
3446        }
3447        nAdd++;
3448        rmatbeg[nAdd] = size;
3449      }
3450      CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3451        rhs, rowSense, rmatbeg,
3452        rmatind, rmatval, NULL);
3453      CPXsetintparam(env, CPX_PARAM_REDUCE,
3454        // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3455        CPX_PREREDUCE_PRIMALONLY);
3456    }
3457    if (getCutoff() < 1.0e50) {
3458      double useCutoff = getCutoff() + offset;
3459      if (bestObjective_ < 1.0e50)
3460        useCutoff = bestObjective_ + offset + 1.0e-7;
3461      cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff * direction);
3462      if (direction > 0.0)
3463        CPXsetdblparam(env, CPX_PARAM_CUTUP, useCutoff); // min
3464      else
3465        CPXsetdblparam(env, CPX_PARAM_CUTLO, useCutoff); // max
3466    }
3467    CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3468    delete[] setSol;
3469    delete[] setVar;
3470    char printBuffer[200];
3471    if (offset) {
3472      sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3473        -offset);
3474      messageHandler()->message(CBC_GENERAL, messages())
3475        << printBuffer << CoinMessageEol;
3476      cpxSolver.setDblParam(OsiObjOffset, offset);
3477    }
3478    cpxSolver.branchAndBound();
3479    double timeTaken = CoinCpuTime() - time1;
3480    sprintf(printBuffer, "Cplex took %g seconds",
3481      timeTaken);
3482    messageHandler()->message(CBC_GENERAL, messages())
3483      << printBuffer << CoinMessageEol;
3484    numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3485    numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3486    double value = cpxSolver.getObjValue() * direction;
3487    if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3488      feasible = true;
3489      clpSolver->setWarmStart(NULL);
3490      // try and do solution
3491      double *newSolution = CoinCopyOfArray(cpxSolver.getColSolution(),
3492        getNumCols());
3493      setBestSolution(CBC_STRONGSOL, value, newSolution);
3494      delete[] newSolution;
3495    }
3496    feasible = false;
3497  }
3498#endif
3499  if (!parentModel_ && (moreSpecialOptions_ & 268435456) != 0) {
3500    // try idiotic idea
3501    CbcObject *obj = new CbcIdiotBranch(this);
3502    obj->setPriority(1); // temp
3503    addObjects(1, &obj);
3504    delete obj;
3505  }
3506
3507  /*
3508      A hook to use clp to quickly explore some part of the tree.
3509    */
3510  if (fastNodeDepth_ == 1000 && /*!parentModel_*/ (specialOptions_ & 2048) == 0) {
3511    fastNodeDepth_ = -1;
3512    CbcObject *obj = new CbcFollowOn(this);
3513    obj->setPriority(1);
3514    addObjects(1, &obj);
3515    delete obj;
3516  }
3517  int saveNumberSolves = numberSolves_;
3518  int saveNumberIterations = numberIterations_;
3519  if ((fastNodeDepth_ >= 0 || (moreSpecialOptions_ & 33554432) != 0)
3520    && /*!parentModel_*/ (specialOptions_ & 2048) == 0) {
3521    // add in a general depth object doClp
3522    int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3523    if ((moreSpecialOptions_ & 33554432) != 0)
3524      type = 12;
3525    else
3526      fastNodeDepth_ += 1000000; // mark as done
3527    CbcObject *obj = new CbcGeneralDepth(this, type);
3528    addObjects(1, &obj);
3529    delete obj;
3530    // fake number of objects
3531    numberObjects_--;
3532    if (parallelMode() < -1) {
3533      // But make sure position is correct
3534      OsiObject *obj2 = object_[numberObjects_];
3535      obj = dynamic_cast<CbcObject *>(obj2);
3536      assert(obj);
3537      obj->setPosition(numberObjects_);
3538    }
3539  }
3540#ifdef COIN_HAS_CLP
3541#ifdef NO_CRUNCH
3542  if (true) {
3543    OsiClpSolverInterface *clpSolver
3544      = dynamic_cast<OsiClpSolverInterface *>(solver_);
3545    if (clpSolver && !parentModel_) {
3546      ClpSimplex *clpSimplex = clpSolver->getModelPtr();
3547      clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
3548      //clpSimplex->startPermanentArrays();
3549      clpSimplex->setPersistenceFlag(2);
3550    }
3551  }
3552#endif
3553#endif
3554  // Save copy of solver
3555  OsiSolverInterface *saveSolver = NULL;
3556  if (!parentModel_ && (specialOptions_ & (512 + 32768)) != 0)
3557    saveSolver = solver_->clone();
3558  double checkCutoffForRestart = 1.0e100;
3559  saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3560  if ((specialOptions_ & 262144) != 0 && !parentModel_) {
3561    // Save stuff and return!
3562    storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3563      solver_->getColLower(),
3564      solver_->getColUpper());
3565    delete[] lowerBefore;
3566    delete[] upperBefore;
3567    delete saveSolver;
3568    if (flipObjective)
3569      flipModel();
3570    return;
3571  }
3572  /*
3573      We've taken the continuous relaxation as far as we can. Time to branch.
3574      The first order of business is to actually create a node. chooseBranch
3575      currently uses strong branching to evaluate branch object candidates,
3576      unless forced back to simple branching. If chooseBranch concludes that a
3577      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3578      == -2) when forced to integer values, it returns here immediately.
3579
3580      Monotone variables trigger a call to resolve(). If the problem remains
3581      feasible, try again to choose a branching variable. At the end of the loop,
3582      resolved == true indicates that some variables were fixed.
3583
3584      Loss of feasibility will result in the deletion of newNode.
3585    */
3586
3587  bool resolved = false;
3588  CbcNode *newNode = NULL;
3589  numberFixedAtRoot_ = 0;
3590  numberFixedNow_ = 0;
3591  if (!parentModel_ && (moreSpecialOptions2_ & 2) != 0) {
3592#ifdef COIN_HAS_CLP
3593    OsiClpSolverInterface *clpSolver
3594      = dynamic_cast<OsiClpSolverInterface *>(solver_);
3595    if (clpSolver) {
3596      if (getCutoff() > 1.0e20) {
3597        printf("Zapping costs\n");
3598        int numberColumns = solver_->getNumCols();
3599        double *zeroCost = new double[numberColumns];
3600        // could make small random
3601        memset(zeroCost, 0, numberColumns * sizeof(double));
3602        solver_->setObjective(zeroCost);
3603        double objValue = solver_->getObjValue();
3604        solver_->setDblParam(OsiObjOffset, -objValue);
3605        clpSolver->getModelPtr()->setObjectiveValue(objValue);
3606        delete[] zeroCost;
3607      } else {
3608        moreSpecialOptions2_ &= ~2;
3609      }
3610    } else {
3611#endif
3612      moreSpecialOptions2_ &= ~2;
3613#ifdef COIN_HAS_CLP
3614    }
3615#endif
3616  }
3617  int numberIterationsAtContinuous = numberIterations_;
3618  //solverCharacteristics_->setSolver(solver_);
3619  if (feasible) {
3620    // mark all cuts as globally valid
3621    int numberCuts = cuts.sizeRowCuts();
3622    resizeWhichGenerator(0, numberCuts);
3623    for (int i = 0; i < numberCuts; i++) {
3624      cuts.rowCutPtr(i)->setGloballyValid();
3625      whichGenerator_[i] = 20000 + (whichGenerator_[i] % 10000);
3626    }
3627#define HOTSTART -1
3628#if HOTSTART < 0
3629    if (bestSolution_ && !parentModel_ && !hotstartSolution_ && (moreSpecialOptions_ & 1024) != 0 && (specialOptions_ & 2048) == 0) {
3630      // Set priorities so only branch on ones we need to
3631      // use djs and see if only few branches needed
3632#ifndef NDEBUG
3633      double integerTolerance = getIntegerTolerance();
3634#endif
3635      bool possible = true;
3636      const double *saveLower = continuousSolver_->getColLower();
3637      const double *saveUpper = continuousSolver_->getColUpper();
3638      for (int i = 0; i < numberObjects_; i++) {
3639        const CbcSimpleInteger *thisOne = dynamic_cast<const CbcSimpleInteger *>(object_[i]);
3640        if (thisOne) {
3641          int iColumn = thisOne->columnNumber();
3642          if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3643            possible = false;
3644            break;
3645          }
3646        } else {
3647          possible = false;
3648          break;
3649        }
3650      }
3651      if (possible) {
3652        OsiSolverInterface *solver = continuousSolver_->clone();
3653        int numberColumns = solver->getNumCols();
3654        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3655          double value = bestSolution_[iColumn];
3656          value = CoinMax(value, saveLower[iColumn]);
3657          value = CoinMin(value, saveUpper[iColumn]);
3658          value = floor(value + 0.5);
3659          if (solver->isInteger(iColumn)) {
3660            solver->setColLower(iColumn, value);
3661            solver->setColUpper(iColumn, value);
3662          }
3663        }
3664        solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3665        // objlim and all slack
3666        double direction = solver->getObjSense();
3667        solver->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
3668        CoinWarmStartBasis *basis = dynamic_cast<CoinWarmStartBasis *>(solver->getEmptyWarmStart());
3669        solver->setWarmStart(basis);
3670        delete basis;
3671        bool changed = true;
3672        hotstartPriorities_ = new int[numberColumns];
3673        for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3674          hotstartPriorities_[iColumn] = 1;
3675        while (changed) {
3676          changed = false;
3677          solver->resolve();
3678          if (!solver->isProvenOptimal()) {
3679            possible = false;
3680            break;
3681          }
3682          const double *dj = solver->getReducedCost();
3683          const double *colLower = solver->getColLower();
3684          const double *colUpper = solver->getColUpper();
3685          const double *solution = solver->getColSolution();
3686          int nAtLbNatural = 0;
3687          int nAtUbNatural = 0;
3688          int nZeroDj = 0;
3689          int nForced = 0;
3690          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3691            double value = solution[iColumn];
3692            value = CoinMax(value, saveLower[iColumn]);
3693            value = CoinMin(value, saveUpper[iColumn]);
3694            if (solver->isInteger(iColumn)) {
3695              assert(fabs(value - solution[iColumn]) <= integerTolerance);
3696              if (hotstartPriorities_[iColumn] == 1) {
3697                if (dj[iColumn] < -1.0e-6) {
3698                  // negative dj
3699                  if (saveUpper[iColumn] == colUpper[iColumn]) {
3700                    nAtUbNatural++;
3701                    hotstartPriorities_[iColumn] = 2;
3702                    solver->setColLower(iColumn, saveLower[iColumn]);
3703                    solver->setColUpper(iColumn, saveUpper[iColumn]);
3704                  } else {
3705                    nForced++;
3706                  }
3707                } else if (dj[iColumn] > 1.0e-6) {
3708                  // positive dj
3709                  if (saveLower[iColumn] == colLower[iColumn]) {
3710                    nAtLbNatural++;
3711                    hotstartPriorities_[iColumn] = 2;
3712                    solver->setColLower(iColumn, saveLower[iColumn]);
3713                    solver->setColUpper(iColumn, saveUpper[iColumn]);
3714                  } else {
3715                    nForced++;
3716                  }
3717                } else {
3718                  // zero dj
3719                  nZeroDj++;
3720                }
3721              }
3722            }
3723          }
3724#if CBC_USEFUL_PRINTING > 1
3725          printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3726            nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3727#endif
3728          if (nAtLbNatural || nAtUbNatural) {
3729            changed = true;
3730          } else {
3731            if (nForced + nZeroDj > 5000 || (nForced + nZeroDj) * 2 > numberIntegers_)
3732              possible = false;
3733          }
3734        }
3735        delete solver;
3736      }
3737      if (possible) {
3738        setHotstartSolution(bestSolution_);
3739        if (!saveCompare) {
3740          // create depth first comparison
3741          saveCompare = nodeCompare_;
3742          // depth first
3743          nodeCompare_ = new CbcCompareDepth();
3744          tree_->setComparison(*nodeCompare_);
3745        }
3746      } else {
3747        delete[] hotstartPriorities_;
3748        hotstartPriorities_ = NULL;
3749      }
3750    }
3751#endif
3752#if HOTSTART > 0
3753    if (hotstartSolution_ && !hotstartPriorities_) {
3754      // Set up hot start
3755      OsiSolverInterface *solver = solver_->clone();
3756      double direction = solver_->getObjSense();
3757      int numberColumns = solver->getNumCols();
3758      double *saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3759      double *saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3760      // move solution
3761      solver->setColSolution(hotstartSolution_);
3762      // point to useful information
3763      const double *saveSolution = testSolution_;
3764      testSolution_ = solver->getColSolution();
3765      OsiBranchingInformation usefulInfo = usefulInformation();
3766      testSolution_ = saveSolution;
3767      /*
3768            Run through the objects and use feasibleRegion() to set variable bounds
3769            so as to fix the variables specified in the objects at their value in this
3770            solution. Since the object list contains (at least) one object for every
3771            integer variable, this has the effect of fixing all integer variables.
3772            */
3773      for (int i = 0; i < numberObjects_; i++)
3774        object_[i]->feasibleRegion(solver, &usefulInfo);
3775      solver->resolve();
3776      assert(solver->isProvenOptimal());
3777      double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0);
3778      const double *dj = solver->getReducedCost();
3779      const double *colLower = solver->getColLower();
3780      const double *colUpper = solver->getColUpper();
3781      const double *solution = solver->getColSolution();
3782      int nAtLbNatural = 0;
3783      int nAtUbNatural = 0;
3784      int nAtLbNaturalZero = 0;
3785      int nAtUbNaturalZero = 0;
3786      int nAtLbFixed = 0;
3787      int nAtUbFixed = 0;
3788      int nAtOther = 0;
3789      int nAtOtherNatural = 0;
3790      int nNotNeeded = 0;
3791      delete[] hotstartSolution_;
3792      hotstartSolution_ = new double[numberColumns];
3793      delete[] hotstartPriorities_;
3794      hotstartPriorities_ = new int[numberColumns];
3795      int *order = (int *)saveUpper;
3796      int nFix = 0;
3797      double bestRatio = COIN_DBL_MAX;
3798      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3799        double value = solution[iColumn];
3800        value = CoinMax(value, saveLower[iColumn]);
3801        value = CoinMin(value, saveUpper[iColumn]);
3802        double sortValue = COIN_DBL_MAX;
3803        if (solver->isInteger(iColumn)) {
3804          assert(fabs(value - solution[iColumn]) <= 1.0e-5);
3805          double value2 = floor(value + 0.5);
3806          if (dj[iColumn] < -1.0e-6) {
3807            // negative dj
3808            //assert (value2==colUpper[iColumn]);
3809            if (saveUpper[iColumn] == colUpper[iColumn]) {
3810              nAtUbNatural++;
3811              sortValue = 0.0;
3812              double value = -dj[iColumn];
3813              if (value > gap)
3814                nFix++;
3815              else if (gap < value * bestRatio)
3816                bestRatio = gap / value;
3817              if (saveLower[iColumn] != colLower[iColumn]) {
3818                nNotNeeded++;
3819                sortValue = 1.0e20;
3820              }
3821            } else if (saveLower[iColumn] == colUpper[iColumn]) {
3822              nAtLbFixed++;
3823              sortValue = dj[iColumn];
3824            } else {
3825              nAtOther++;
3826              sortValue = 0.0;
3827              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
3828                nNotNeeded++;
3829                sortValue = 1.0e20;
3830              }
3831            }
3832          } else if (dj[iColumn] > 1.0e-6) {
3833            // positive dj
3834            //assert (value2==colLower[iColumn]);
3835            if (saveLower[iColumn] == colLower[iColumn]) {
3836              nAtLbNatural++;
3837              sortValue = 0.0;
3838              double value = dj[iColumn];
3839              if (value > gap)
3840                nFix++;
3841              else if (gap < value * bestRatio)
3842                bestRatio = gap / value;
3843              if (saveUpper[iColumn] != colUpper[iColumn]) {
3844                nNotNeeded++;
3845                sortValue = 1.0e20;
3846              }
3847            } else if (saveUpper[iColumn] == colLower[iColumn]) {
3848              nAtUbFixed++;
3849              sortValue = -dj[iColumn];
3850            } else {
3851              nAtOther++;
3852              sortValue = 0.0;
3853              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
3854                nNotNeeded++;
3855                sortValue = 1.0e20;
3856              }
3857            }
3858          } else {
3859            // zero dj
3860            if (value2 == saveUpper[iColumn]) {
3861              nAtUbNaturalZero++;
3862              sortValue = 0.0;
3863              if (saveLower[iColumn] != colLower[iColumn]) {
3864                nNotNeeded++;
3865                sortValue = 1.0e20;
3866              }
3867            } else if (value2 == saveLower[iColumn]) {
3868              nAtLbNaturalZero++;
3869              sortValue = 0.0;
3870            } else {
3871              nAtOtherNatural++;
3872              sortValue = 0.0;
3873              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
3874                nNotNeeded++;
3875                sortValue = 1.0e20;
3876              }
3877            }
3878          }
3879#if HOTSTART == 3
3880          sortValue = -fabs(dj[iColumn]);
3881#endif
3882        }
3883        hotstartSolution_[iColumn] = value;
3884        saveLower[iColumn] = sortValue;
3885        order[iColumn] = iColumn;
3886      }
3887      COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3888        nFix, bestRatio, gap));
3889      int nNeg = 0;
3890      CoinSort_2(saveLower, saveLower + numberColumns, order);
3891      for (int i = 0; i < numberColumns; i++) {
3892        if (saveLower[i] < 0.0) {
3893          nNeg++;
3894#if HOTSTART == 2 || HOTSTART == 3
3895          // swap sign ?
3896          saveLower[i] = -saveLower[i];
3897#endif
3898        }
3899      }
3900      CoinSort_2(saveLower, saveLower + nNeg, order);
3901      for (int i = 0; i < numberColumns; i++) {
3902#if HOTSTART == 1
3903        hotstartPriorities_[order[i]] = 100;
3904#else
3905          hotstartPriorities_[order[i]] = -(i + 1);
3906#endif
3907      }
3908      COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3909        nAtLbNatural,
3910        nAtUbNatural,
3911        nAtLbNaturalZero,
3912        nAtUbNaturalZero,
3913        nAtLbFixed,
3914        nAtUbFixed,
3915        nAtOther,
3916        nAtOtherNatural, nNotNeeded, nNeg));
3917      delete[] saveLower;
3918      delete[] saveUpper;
3919      if (!saveCompare) {
3920        // create depth first comparison
3921        saveCompare = nodeCompare_;
3922        // depth first
3923        nodeCompare_ = new CbcCompareDepth();
3924        tree_->setComparison(*nodeCompare_);
3925      }
3926    }
3927#endif
3928    newNode = new CbcNode;
3929    // Set objective value (not so obvious if NLP etc)
3930    setObjectiveValue(newNode, NULL);
3931    anyAction = -1;
3932    // To make depth available we may need a fake node
3933    CbcNode fakeNode;
3934    if (!currentNode_) {
3935      // Not true if sub trees assert (!numberNodes_);
3936      currentNode_ = &fakeNode;
3937    }
3938    phase_ = 3;
3939    // only allow 1000 passes
3940    int numberPassesLeft = 1000;
3941    // This is first crude step
3942    if (numberAnalyzeIterations_ && !parentModel_) {
3943      delete[] analyzeResults_;
3944      //int numberColumns = solver_->getNumCols();
3945      analyzeResults_ = new double[5 * numberIntegers_];
3946      numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3947      if (numberFixedAtRoot_ > 0) {
3948        COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3949        setPointers(solver_);
3950        numberFixedNow_ = numberFixedAtRoot_;
3951      } else if (numberFixedAtRoot_ < 0) {
3952        COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3953        anyAction = -2;
3954        delete newNode;
3955        newNode = NULL;
3956        feasible = false;
3957      }
3958    }
3959    OsiSolverBranch *branches = NULL;
3960    if (feasible)
3961      anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3962        NULL, NULL, NULL, branches);
3963    if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3964      if (anyAction != -2) {
3965        // zap parent nodeInfo
3966#ifdef COIN_DEVELOP
3967        printf("zapping CbcNodeInfo %x\n", newNode->nodeInfo()->parent());
3968#endif
3969        if (newNode->nodeInfo())
3970          newNode->nodeInfo()->nullParent();
3971      }
3972      delete newNode;
3973      newNode = NULL;
3974      feasible = false;
3975    }
3976  }
3977  if (newNode && probingInfo_) {
3978    int number01 = probingInfo_->numberIntegers();
3979    //const fixEntry * entry = probingInfo_->fixEntries();
3980    const int *toZero = probingInfo_->toZero();
3981    //const int * toOne = probingInfo_->toOne();
3982    //const int * integerVariable = probingInfo_->integerVariable();
3983    if (toZero[number01]) {
3984      CglTreeProbingInfo info(*probingInfo_);
3985      if ((moreSpecialOptions2_ & 64) != 0 && !parentModel_) {
3986        /*
3987                Marginal idea. Further exploration probably good. Build some extra
3988                cliques from probing info. Not quite worth the effort?
3989              */
3990        CglProbing generator1;
3991        generator1.setUsingObjective(false);
3992        generator1.setMaxPass(1);
3993        generator1.setMaxPassRoot(1);
3994        generator1.setMaxLook(100);
3995        generator1.setRowCuts(3);
3996        generator1.setMaxElements(300);
3997        generator1.setMaxProbeRoot(solver_->getNumCols());
3998        CoinThreadRandom randomGenerator;
3999        //CglTreeProbingInfo info(solver_);
4000        info.level = 0;
4001        info.formulation_rows = solver_->getNumRows();
4002        info.inTree = false;
4003        if (parentModel_) {
4004          info.parentSolver = parentModel_->continuousSolver();
4005          // indicate if doing full search
4006          info.hasParent = ((specialOptions_ & 67108864) == 0) ? 1 : 2;
4007        } else {
4008          info.hasParent = 0;
4009          info.parentSolver = NULL;
4010        }
4011        info.originalColumns = originalColumns();
4012        info.randomNumberGenerator = &randomGenerator;
4013        info.pass = 4;
4014        generator1.setMode(8);
4015        OsiCuts cs;
4016        generator1.generateCutsAndModify(*solver_, cs, &info);
4017        // very clunky
4018        OsiSolverInterface *temp = generator1.cliqueModel(solver_, 2);
4019        CglPreProcess dummy;
4020        OsiSolverInterface *newSolver = dummy.cliqueIt(*temp, 0.0001);
4021        delete temp;
4022        OsiSolverInterface *fake = NULL;
4023        if (newSolver) {
4024#if 0
4025                int numberCliques = generator1.numberCliques();
4026                cliqueEntry * entry = generator1.cliqueEntry();
4027                cliqueType * type = new cliqueType [numberCliques];
4028                int * start = new int [numberCliques+1];
4029                start[numberCliques]=2*numberCliques;
4030                int n=0;
4031                for (int i=0;i<numberCliques;i++) {
4032                  start[i]=2*i;
4033                  setOneFixesInCliqueEntry(entry[2*i],true);
4034                  setOneFixesInCliqueEntry(entry[2*i+1],true);
4035                  type[i]=0;
4036                }
4037                fake = info.analyze(*solver_, 1,numberCliques,start,
4038                                    entry,type);
4039                delete [] type;
4040                delete [] entry;
4041#else
4042        fake = info.analyze(*newSolver, 1, -1);
4043#endif
4044          delete newSolver;
4045        } else {
4046          fake = info.analyze(*solver_, 1);
4047        }
4048        if (fake) {
4049          //fake->writeMps("fake");
4050          CglFakeClique cliqueGen(fake);
4051          cliqueGen.setStarCliqueReport(false);
4052          cliqueGen.setRowCliqueReport(false);
4053          cliqueGen.setMinViolation(0.1);
4054          addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
4055          generator_[numberCutGenerators_ - 1]->setTiming(true);
4056          for (int i = 0; i < numberCutGenerators_; i++) {
4057            CglKnapsackCover *cutGen = dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
4058            if (cutGen) {
4059              cutGen->createCliques(*fake, 2, 200, false);
4060            }
4061          }
4062        }
4063      }
4064      if (probingInfo_->packDown()) {
4065#if CBC_USEFUL_PRINTING > 1
4066        printf("%d implications on %d 0-1\n", toZero[number01], number01);
4067#endif
4068        // Create a cut generator that remembers implications discovered at root.
4069        CglImplication implication(probingInfo_);
4070        addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
4071        generator_[numberCutGenerators_ - 1]->setGlobalCuts(true);
4072        generator_[numberCutGenerators_ - 1]->setTiming(true);
4073      } else {
4074        delete probingInfo_;
4075        probingInfo_ = NULL;
4076      }
4077    } else {
4078      delete probingInfo_;
4079
4080      probingInfo_ = NULL;
4081    }
4082  }
4083  /*
4084      At this point, the root subproblem is infeasible or fathomed by bound
4085      (newNode == NULL), or we're live with an objective value that satisfies the
4086      current objective cutoff.
4087    */
4088  assert(!newNode || newNode->objectiveValue() <= cutoff);
4089  // Save address of root node as we don't want to delete it
4090  /*
4091      The common case is that the lp relaxation is feasible but doesn't satisfy
4092      integrality (i.e., newNode->branchingObject(), indicating we've been able to
4093      select a branching variable). Remove any cuts that have gone slack due to
4094      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
4095      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
4096      addCuts()). If, by some miracle, we have an integral solution at the root
4097      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
4098      a valid solution for use by setBestSolution().
4099    */
4100  CoinWarmStartBasis *lastws = NULL;
4101  if (feasible && newNode->branchingObject()) {
4102    if (resolved) {
4103      takeOffCuts(cuts, false, NULL);
4104#ifdef CHECK_CUT_COUNTS
4105      {
4106        printf("Number of rows after chooseBranch fix (root)"
4107               "(active only) %d\n",
4108          numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_);
4109        const CoinWarmStartBasis *debugws = dynamic_cast<const CoinWarmStartBasis *>(solver_->getWarmStart());
4110        debugws->print();
4111        delete debugws;
4112      }
4113#endif
4114    }
4115    //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
4116    //newNode->nodeInfo()->addCuts(cuts,
4117    //                   newNode->numberBranches(),whichGenerator_) ;
4118    if (lastws)
4119      delete lastws;
4120    lastws = dynamic_cast<CoinWarmStartBasis *>(solver_->getWarmStart());
4121  }
4122  /*
4123      Continuous data to be used later
4124    */
4125  continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4126  continuousInfeasibilities_ = 0;
4127  if (newNode) {
4128    continuousObjective_ = newNode->objectiveValue();
4129    delete[] continuousSolution_;
4130    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4131      numberColumns);
4132    continuousInfeasibilities_ = newNode->numberUnsatisfied();
4133  }
4134  /*
4135      Bound may have changed so reset in objects
4136    */
4137  {
4138    int i;
4139    for (i = 0; i < numberObjects_; i++)
4140      object_[i]->resetBounds(solver_);
4141  }
4142  /*
4143      Feasible? Then we should have either a live node prepped for future
4144      expansion (indicated by variable() >= 0), or (miracle of miracles) an
4145      integral solution at the root node.
4146
4147      initializeInfo sets the reference counts in the nodeInfo object.  Since
4148      this node is still live, push it onto the heap that holds the live set.
4149    */
4150  if (newNode) {
4151    if (newNode->branchingObject()) {
4152      newNode->initializeInfo();
4153      tree_->push(newNode);
4154      // save pointer to root node - so can pick up bounds
4155      if (!topOfTree_)
4156        topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo());
4157      if (statistics_) {
4158        if (numberNodes2_ == maximumStatistics_) {
4159          maximumStatistics_ = 2 * maximumStatistics_;
4160          CbcStatistics **temp = new CbcStatistics *[maximumStatistics_];
4161          memset(temp, 0, maximumStatistics_ * sizeof(CbcStatistics *));
4162          memcpy(temp, statistics_, numberNodes2_ * sizeof(CbcStatistics *));
4163          delete[] statistics_;
4164          statistics_ = temp;
4165        }
4166        assert(!statistics_[numberNodes2_]);
4167        statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4168      }
4169      numberNodes2_++;
4170#ifdef CHECK_NODE
4171      printf("Node %x on tree\n", newNode);
4172#endif
4173    } else {
4174      // continuous is integer
4175      double objectiveValue = newNode->objectiveValue();
4176      setBestSolution(CBC_SOLUTION, objectiveValue,
4177        solver_->getColSolution());
4178      if (eventHandler) {
4179        // we are stopping anyway so no need to test return code
4180        eventHandler->event(CbcEventHandler::solution);
4181      }
4182      delete newNode;
4183      newNode = NULL;
4184    }
4185  }
4186
4187  if (printFrequency_ <= 0) {
4188    printFrequency_ = 1000;
4189    if (getNumCols() > 2000)
4190      printFrequency_ = 100;
4191  }
4192  /*
4193      It is possible that strong branching fixes one variable and then the code goes round
4194      again and again.  This can take too long.  So we need to warn user - just once.
4195    */
4196  numberLongStrong_ = 0;
4197  CbcNode *createdNode = NULL;
4198#ifdef CBC_THREAD
4199  if ((specialOptions_ & 2048) != 0)
4200    numberThreads_ = 0;
4201  if (numberThreads_) {
4202    nodeCompare_->sayThreaded(); // need to use addresses
4203    master_ = new CbcBaseModel(*this,
4204      (parallelMode() < -1) ? 1 : 0);
4205    masterThread_ = master_->masterThread();
4206  }
4207#endif
4208#ifdef COIN_HAS_CLP
4209  {
4210    OsiClpSolverInterface *clpSolver
4211      = dynamic_cast<OsiClpSolverInterface *>(solver_);
4212    if (clpSolver && !parentModel_) {
4213      clpSolver->computeLargestAway();
4214    }
4215  }
4216#endif
4217  /*
4218      At last, the actual branch-and-cut search loop, which will iterate until
4219      the live set is empty or we hit some limit (integrality gap, time, node
4220      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4221      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4222      add the node to the live set.
4223
4224      The first action is to winnow the live set to remove nodes which are worse
4225      than the current objective cutoff.
4226    */
4227  if (solver_->getRowCutDebuggerAlways()) {
4228    OsiRowCutDebugger *debuggerX = const_cast<OsiRowCutDebugger *>(solver_->getRowCutDebuggerAlways());
4229    const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger();
4230    if (!debugger) {
4231      // infeasible!!
4232      printf("before search\n");
4233      const double *lower = solver_->getColLower();
4234      const double *upper = solver_->getColUpper();
4235      const double *solution = debuggerX->optimalSolution();
4236      int numberColumns = solver_->getNumCols();
4237      for (int i = 0; i < numberColumns; i++) {
4238        if (solver_->isInteger(i)) {
4239          if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
4240            printf("**** ");
4241          printf("%d %g <= %g <= %g\n",
4242            i, lower[i], solution[i], upper[i]);
4243        }
4244      }
4245      //abort();
4246    }
4247  }
4248  {
4249    // may be able to change cutoff now
4250    double cutoff = getCutoff();
4251    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
4252    if (cutoff > bestObjective_ - increment) {
4253      cutoff = bestObjective_ - increment;
4254      setCutoff(cutoff);
4255    }
4256  }
4257#ifdef CBC_THREAD
4258  bool goneParallel = false;
4259#endif
4260#define MAX_DEL_NODE 1
4261  CbcNode *delNode[MAX_DEL_NODE + 1];
4262  int nDeleteNode = 0;
4263  // For Printing etc when parallel
4264  int lastEvery1000 = 0;
4265  int lastPrintEvery = 0;
4266  int numberConsecutiveInfeasible = 0;
4267#define PERTURB_IN_FATHOM
4268#ifdef PERTURB_IN_FATHOM
4269  // allow in fathom
4270  if ((moreSpecialOptions_ & 262144) != 0)
4271    specialOptions_ |= 131072;
4272#endif
4273  while (true) {
4274    lockThread();
4275#ifdef COIN_HAS_CLP
4276    // See if we want dantzig row choice
4277    goToDantzig(100, savePivotMethod);
4278#endif
4279    //#define REPORT_DYNAMIC 2
4280#if REPORT_DYNAMIC
4281    if (numberNodes_ && !parentModel_ && (tree_->empty() || (numberNodes_ % 10000) == 0)) {
4282      // Get average up and down costs
4283      double averageUp = 0.0;
4284      double averageDown = 0.0;
4285      int numberUp = 0;
4286      int numberDown = 0;
4287      int minTimesDown = COIN_INT_MAX;
4288      int maxTimesDown = 0;
4289      int neverBranchedDown = 0;
4290      int infeasibleTimesDown = 0;
4291      int minTimesUp = COIN_INT_MAX;
4292      int maxTimesUp = 0;
4293      int infeasibleTimesUp = 0;
4294      int neverBranchedUp = 0;
4295      int neverBranched = 0;
4296      int i;
4297      int numberInts = 0;
4298      bool endOfSearch = tree_->empty();
4299      int numberUp2 = 0;
4300      int numberDown2 = 0;
4301      for (i = 0; i < numberObjects_; i++) {
4302        OsiObject *object = object_[i];
4303        CbcSimpleIntegerDynamicPseudoCost *dynamicObject = dynamic_cast<CbcSimpleIntegerDynamicPseudoCost *>(object);
4304        if (dynamicObject) {
4305          numberInts++;
4306          if (dynamicObject->numberTimesUp() || dynamicObject->numberTimesDown()) {
4307            int nUp = 0;
4308            int nDown = 0;
4309            double up = 0.0;
4310            double down = 0.0;
4311            if (dynamicObject->numberTimesUp()) {
4312              numberUp++;
4313              nUp = dynamicObject->numberTimesUp();
4314              minTimesUp = CoinMin(minTimesUp, nUp);
4315              maxTimesUp = CoinMax(maxTimesUp, nUp);
4316              up = dynamicObject->upDynamicPseudoCost();
4317              averageUp += up;
4318              numberUp2 += nUp;
4319              infeasibleTimesUp += dynamicObject->numberTimesUpInfeasible();
4320            } else {
4321              neverBranchedUp++;
4322            }
4323            if (dynamicObject->numberTimesDown()) {
4324              numberDown++;
4325              nDown = dynamicObject->numberTimesDown();
4326              minTimesDown = CoinMin(minTimesDown, nDown);
4327              maxTimesDown = CoinMax(maxTimesDown, nDown);
4328              down = dynamicObject->downDynamicPseudoCost();
4329              averageDown += down;
4330              numberDown2 += dynamicObject->numberTimesDown();
4331              infeasibleTimesDown += dynamicObject->numberTimesDownInfeasible();
4332            } else {
4333              neverBranchedDown++;
4334            }
4335#if REPORT_DYNAMIC > 1
4336#if REPORT_DYNAMIC == 2
4337            if (endOfSearch && numberIntegers_ < 400) {
4338#elif REPORT_DYNAMIC == 3
4339            if (endOfSearch) {
4340#else
4341            {
4342#endif
4343              dynamicObject->print(0, 0.0);
4344            }
4345#endif
4346          } else {
4347            neverBranched++;
4348#if REPORT_DYNAMIC > 2
4349#if REPORT_DYNAMIC == 3
4350            if (endOfSearch && numberIntegers_ < 400) {
4351#elif REPORT_DYNAMIC == 4
4352            if (endOfSearch) {
4353#else
4354            {
4355#endif
4356              printf("col %d - never branched on\n", dynamicObject->columnNumber());
4357            }
4358#endif
4359          }
4360        }
4361      }
4362      if (numberUp)
4363        averageUp /= static_cast<double>(numberUp);
4364      else
4365        averageUp = 0.0;
4366      if (numberDown)
4367        averageDown /= static_cast<double>(numberDown);
4368      else
4369        averageDown = 0.0;
4370      printf("Report for %d variables (%d never branched on) after %d nodes - total solves down %d up %d\n",
4371        numberInts, neverBranched, numberNodes_, numberDown2, numberUp2);
4372      if ((neverBranchedDown || neverBranchedUp) && endOfSearch)
4373        printf("odd %d never branched down and %d never branched up\n",
4374          neverBranchedDown, neverBranchedUp);
4375      printf("down average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
4376        static_cast<double>(numberDown2) / numberDown, infeasibleTimesDown, averageDown,
4377        minTimesDown, maxTimesDown);
4378      printf("up average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
4379        static_cast<double>(numberUp2) / numberUp, infeasibleTimesUp, averageUp,
4380        minTimesUp, maxTimesUp);
4381    }
4382#endif
4383    if (tree_->empty()) {
4384#ifdef CBC_THREAD
4385      if (parallelMode() > 0 && master_) {
4386        int anyLeft = master_->waitForThreadsInTree(0);
4387        if (!anyLeft) {
4388          master_->stopThreads(-1);
4389          break;
4390        }
4391      } else {
4392        break;
4393      }
4394#else
4395    break;
4396#endif
4397    } else {
4398      unlockThread();
4399    }
4400    // If done 50/100 nodes see if worth trying reduction
4401    if (numberNodes_ >= nextCheckRestart) {
4402      if (nextCheckRestart < 100)
4403        nextCheckRestart = 100;
4404      else
4405        nextCheckRestart = COIN_INT_MAX;
4406#ifdef COIN_HAS_CLP
4407      OsiClpSolverInterface *clpSolver
4408        = dynamic_cast<OsiClpSolverInterface *>(solver_);
4409      if (clpSolver && ((specialOptions_ & 131072) == 0) && true) {
4410        ClpSimplex *simplex = clpSolver->getModelPtr();
4411        int perturbation = simplex->perturbation();
4412#if CBC_USEFUL_PRINTING > 1
4413        printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4414          numberIterations_, saveNumberIterations,
4415          numberSolves_, saveNumberSolves, perturbation);
4416#endif
4417        if (perturbation == 50 && (numberIterations_ - saveNumberIterations) < 8 * (numberSolves_ - saveNumberSolves)) {
4418          // switch off perturbation
4419          simplex->setPerturbation(100);
4420#if CBC_USEFUL_PRINTING > 1
4421          printf("Perturbation switched off\n");
4422#endif
4423        }
4424      }
4425#endif
4426      /*
4427              Decide if we want to do a restart.
4428            */
4429      if (saveSolver && (specialOptions_ & (512 + 32768)) != 0) {
4430        bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() && (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4431        int numberColumns = getNumCols();
4432        if (tryNewSearch) {
4433          // adding increment back allows current best - tiny bit weaker
4434          checkCutoffForRestart = getCutoff() + getCutoffIncrement();
4435#if CBC_USEFUL_PRINTING > 1
4436          printf("after %d nodes, cutoff %g - looking\n",
4437            numberNodes_, getCutoff());
4438#endif
4439          saveSolver->resolve();
4440          double direction = saveSolver->getObjSense();
4441          double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction;
4442          double tolerance;
4443          saveSolver->getDblParam(OsiDualTolerance, tolerance);
4444          if (gap <= 0.0)
4445            gap = tolerance;
4446          gap += 100.0 * tolerance;
4447          double integerTolerance = getDblParam(CbcIntegerTolerance);
4448
4449          const double *lower = saveSolver->getColLower();
4450          const double *upper = saveSolver->getColUpper();
4451          const double *solution = saveSolver->getColSolution();
4452          const double *reducedCost = saveSolver->getReducedCost();
4453
4454          int numberFixed = 0;
4455          int numberFixed2 = 0;
4456#ifdef COIN_DEVELOP
4457          printf("gap %g\n", gap);
4458#endif
4459          for (int i = 0; i < numberIntegers_; i++) {
4460            int iColumn = integerVariable_[i];
4461            double djValue = direction * reducedCost[iColumn];
4462            if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4463              if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4464                //printf("%d to lb on dj of %g - bounds %g %g\n",
4465                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4466                saveSolver->setColUpper(iColumn, lower[iColumn]);
4467                numberFixed++;
4468              } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4469                //printf("%d to ub on dj of %g - bounds %g %g\n",
4470                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4471                saveSolver->setColLower(iColumn, upper[iColumn]);
4472                numberFixed++;
4473              }
4474            } else {
4475              //printf("%d has dj of %g - already fixed to %g\n",
4476              //     iColumn,djValue,lower[iColumn]);
4477              numberFixed2++;
4478            }
4479          }
4480#ifdef COIN_DEVELOP
4481          if ((specialOptions_ & 1) != 0) {
4482            const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger();
4483            if (debugger) {
4484              printf("Contains optimal\n");
4485              OsiSolverInterface *temp = saveSolver->clone();
4486              const double *solution = debugger->optimalSolution();
4487              const double *lower = temp->getColLower();
4488              const double *upper = temp->getColUpper();
4489              int n = temp->getNumCols();
4490              for (int i = 0; i < n; i++) {
4491                if (temp->isInteger(i)) {
4492                  double value = floor(solution[i] + 0.5);
4493                  assert(value >= lower[i] && value <= upper[i]);
4494                  temp->setColLower(i, value);
4495                  temp->setColUpper(i, value);
4496                }
4497              }
4498              temp->writeMps("reduced_fix");
4499              delete temp;
4500              saveSolver->writeMps("reduced");
4501            } else {
4502              abort();
4503            }
4504          }
4505          printf("Restart could fix %d integers (%d already fixed)\n",
4506            numberFixed + numberFixed2, numberFixed2);
4507#endif
4508          numberFixed += numberFixed2;
4509          if (numberFixed * 10 < numberColumns && numberFixed * 4 < numberIntegers_)
4510            tryNewSearch = false;
4511        }
4512#ifdef CONFLICT_CUTS
4513        // temporary
4514        //if ((moreSpecialOptions_&4194304)!=0)
4515        //tryNewSearch=false;
4516#endif
4517        if (tryNewSearch) {
4518          // back to solver without cuts?
4519          OsiSolverInterface *solver2 = saveSolver->clone();
4520          const double *lower = saveSolver->getColLower();
4521          const double *upper = saveSolver->getColUpper();
4522          for (int i = 0; i < numberIntegers_; i++) {
4523            int iColumn = integerVariable_[i];
4524            solver2->setColLower(iColumn, lower[iColumn]);
4525            solver2->setColUpper(iColumn, upper[iColumn]);
4526          }
4527          // swap
4528          delete saveSolver;
4529          saveSolver = solver2;
4530          double *newSolution = new double[numberColumns];
4531          double objectiveValue = checkCutoffForRestart;
4532          // Save the best solution so far.
4533          CbcSerendipity heuristic(*this);
4534          if (bestSolution_)
4535            heuristic.setInputSolution(bestSolution_, bestObjective_);
4536          // Magic number
4537          heuristic.setFractionSmall(0.8);
4538          // `pumpTune' to stand-alone solver for explanations.
4539          heuristic.setFeasibilityPumpOptions(1008013);
4540          // Use numberNodes to say how many are original rows
4541          heuristic.setNumberNodes(continuousSolver_->getNumRows());
4542#ifdef COIN_DEVELOP
4543          if (continuousSolver_->getNumRows() < solver_->getNumRows())
4544            printf("%d rows added ZZZZZ\n",
4545              solver_->getNumRows() - continuousSolver_->getNumRows());
4546#endif
4547          int returnCode = heuristic.smallBranchAndBound(saveSolver,
4548            -1, newSolution,
4549            objectiveValue,
4550            checkCutoffForRestart, "Reduce");
4551          if (returnCode < 0) {
4552#ifdef COIN_DEVELOP
4553            printf("Restart - not small enough to do search after fixing\n");
4554#endif
4555            delete[] newSolution;
4556          } else {
4557            // 1 for sol'n, 2 for finished, 3 for both
4558            if ((returnCode & 1) != 0) {
4559              // increment number of solutions so other heuristics can test
4560              numberSolutions_++;
4561              numberHeuristicSolutions_++;
4562              lastHeuristic_ = NULL;
4563              setBestSolution(CBC_ROUNDING, objectiveValue, newSolution);
4564            }
4565            delete[] newSolution;
4566#ifdef CBC_THREAD
4567            if (master_) {
4568              lockThread();
4569              if (parallelMode() > 0) {
4570                while (master_->waitForThreadsInTree(0)) {
4571                  lockThread();
4572                  double dummyBest;
4573                  tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4574                  //unlockThread();
4575                }
4576              } else {
4577                double dummyBest;
4578                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4579              }
4580              master_->waitForThreadsInTree(2);
4581              delete master_;
4582              master_ = NULL;
4583              masterThread_ = NULL;
4584            }
4585#endif
4586            if (tree_->size()) {
4587              double dummyBest;
4588              tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4589            }
4590            break;
4591          }
4592        }
4593        delete saveSolver;
4594        saveSolver = NULL;
4595      }
4596    }
4597    /*
4598          Check for abort on limits: node count, solution count, time, integrality gap.
4599        */
4600    if (!(numberNodes_ < intParam_[CbcMaxNumNode] && numberSolutions_ < intParam_[CbcMaxNumSol] && !maximumSecondsReached() && !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))) {
4601      // out of loop
4602      break;
4603    }
4604#ifdef BONMIN
4605    assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4606#endif
4607// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4608// we need to replace the comparison function to dive, and that requires reordering to retain the
4609// heap property.
4610#define DIVE_WHEN 1000
4611#define DIVE_STOP 2000
4612    int kNode = numberNodes_ % 4000;
4613    if (numberNodes_ < 100000 && kNode > DIVE_WHEN && kNode <= DIVE_STOP) {
4614      if (!parallelMode()) {
4615        if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4616          CbcCompareDefault *compare = dynamic_cast<CbcCompareDefault *>(nodeCompare_);
4617          // Don't interfere if user has replaced the compare function.
4618          if (compare) {
4619            //printf("Redoing tree\n");
4620            compare->startDive(this);
4621            numberConsecutiveInfeasible = 0;
4622          }
4623        }
4624      }
4625    }
4626    // replace current cutoff?
4627    if (cutoff > getCutoff()) {
4628      double newCutoff = getCutoff();
4629      if (analyzeResults_) {
4630        // see if we could fix any (more)
4631        int n = 0;
4632        double *newLower = analyzeResults_;
4633        double *objLower = newLower + numberIntegers_;
4634        double *newUpper = objLower + numberIntegers_;
4635        double *objUpper = newUpper + numberIntegers_;
4636        for (int i = 0; i < numberIntegers_; i++) {
4637          if (objLower[i] > newCutoff) {
4638            n++;
4639            if (objUpper[i] > newCutoff) {
4640              newCutoff = -COIN_DBL_MAX;
4641              break;
4642            }
4643            // add as global cut
4644            objLower[i] = -COIN_DBL_MAX;
4645            OsiRowCut rc;
4646            rc.setLb(newLower[i]);
4647            rc.setUb(COIN_DBL_MAX);
4648            double one = 1.0;
4649            rc.setRow(1, integerVariable_ + i, &one, false);
4650            rc.setGloballyValidAsInteger(2);
4651            globalCuts_.addCutIfNotDuplicate(rc);
4652          } else if (objUpper[i] > newCutoff) {
4653            n++;
4654            // add as global cut
4655            objUpper[i] = -COIN_DBL_MAX;
4656            OsiRowCut rc;
4657            rc.setLb(-COIN_DBL_MAX);
4658            rc.setUb(newUpper[i]);
4659            double one = 1.0;
4660            rc.setRow(1, integerVariable_ + i, &one, false);
4661            rc.setGloballyValidAsInteger(2);
4662            globalCuts_.addCutIfNotDuplicate(rc);
4663          }
4664        }
4665        if (newCutoff == -COIN_DBL_MAX) {
4666          COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4667        } else if (n > numberFixedNow_) {
4668          COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4669          numberFixedNow_ = n;
4670        }
4671      }
4672      if (eventHandler) {
4673        if (!eventHandler->event(CbcEventHandler::solution)) {
4674          eventHappened_ = true; // exit
4675        }
4676        newCutoff = getCutoff();
4677      }
4678      lockThread();
4679      /*
4680              Clean the tree to reflect the new solution, then see if the
4681              node comparison predicate wants to make any changes. If so,
4682              call setComparison for the side effect of rebuilding the heap.
4683            */
4684      tree_->cleanTree(this, newCutoff, bestPossibleObjective_);
4685      if (nodeCompare_->newSolution(this) || nodeCompare_->newSolution(this, continuousObjective_, continuousInfeasibilities_)) {
4686        tree_->setComparison(*nodeCompare_);
4687      }
4688      if (tree_->empty()) {
4689        continue;
4690      }
4691      unlockThread();
4692    }
4693    cutoff = getCutoff();
4694    /*
4695            Periodic activities: Opportunities to
4696            + tweak the nodeCompare criteria,
4697            + check if we've closed the integrality gap enough to quit,
4698            + print a summary line to let the user know we're working
4699        */
4700    if (numberNodes_ >= lastEvery1000) {
4701      lockThread();
4702#ifdef COIN_HAS_CLP
4703      // See if we want dantzig row choice
4704      goToDantzig(1000, savePivotMethod);
4705#endif
4706      lastEvery1000 = numberNodes_ + 1000;
4707      bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_);
4708#ifdef CHECK_CUT_SIZE
4709      verifyCutSize(tree_, *this);
4710#endif
4711      // redo tree if requested
4712      if (redoTree)
4713        tree_->setComparison(*nodeCompare_);
4714      unlockThread();
4715    }
4716    // Had hotstart before, now switched off
4717    if (saveCompare && !hotstartSolution_) {
4718      // hotstart switched off
4719      delete nodeCompare_; // off depth first
4720      nodeCompare_ = saveCompare;
4721      saveCompare = NULL;
4722      // redo tree
4723      lockThread();
4724      tree_->setComparison(*nodeCompare_);
4725      unlockThread();
4726    }
4727    if (numberNodes_ >= lastPrintEvery) {
4728      lastPrintEvery = numberNodes_ + printFrequency_;
4729      lockThread();
4730      int nNodes = tree_->size();
4731
4732      //MODIF PIERRE
4733      bestPossibleObjective_ = tree_->getBestPossibleObjective();
4734#ifdef CBC_THREAD
4735      if (parallelMode() > 0 && master_) {
4736        // need to adjust for ones not on tree
4737        int numberThreads = master_->numberThreads();
4738        for (int i = 0; i < numberThreads; i++) {
4739          CbcThread *child = master_->child(i);
4740          if (child->node()) {
4741            // adjust
4742            double value = child->node()->objectiveValue();
4743            bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
4744          }
4745        }
4746      }
4747#endif
4748      unlockThread();
4749#if CBC_USEFUL_PRINTING > 1
4750      if (getCutoff() < 1.0e20) {
4751        if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 && !parentModel_)
4752          printf("model cutoff in status %g, best %g, increment %g\n",
4753            getCutoff(), bestObjective_, getCutoffIncrement());
4754        assert(getCutoff() < bestObjective_ - getCutoffIncrement() + 1.0e-6 + 1.0e-10 * fabs(bestObjective_));
4755      }
4756#endif
4757      if (!intParam_[CbcPrinting]) {
4758        // Parallel may not have any nodes
4759        if (!nNodes)
4760          bestPossibleObjective_ = lastBestPossibleObjective;
4761        else
4762          lastBestPossibleObjective = bestPossibleObjective_;
4763        messageHandler()->message(CBC_STATUS, messages())
4764          << numberNodes_ << CoinMax(nNodes, 1) << bestObjective_ << bestPossibleObjective_
4765          << getCurrentSeconds()
4766          << CoinMessageEol;
4767      } else if (intParam_[CbcPrinting] == 1) {
4768        messageHandler()->message(CBC_STATUS2, messages())
4769          << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4770          << tree_->lastDepth() << tree_->lastUnsatisfied()
4771          << tree_->lastObjective() << numberIterations_
4772          << getCurrentSeconds()
4773          << CoinMessageEol;
4774      } else if (!numberExtraIterations_) {
4775        messageHandler()->message(CBC_STATUS2, messages())
4776          << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4777          << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4778          << getCurrentSeconds()
4779          << CoinMessageEol;
4780      } else {
4781        messageHandler()->message(CBC_STATUS3, messages())
4782          << numberNodes_ << numberFathoms_ << numberExtraNodes_ << nNodes
4783          << bestObjective_ << bestPossibleObjective_
4784          << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4785          << getCurrentSeconds()
4786          << CoinMessageEol;
4787      }
4788#ifdef COIN_HAS_NTY
4789      if (symmetryInfo_)
4790        symmetryInfo_->statsOrbits(this, 1);
4791#endif
4792#if PRINT_CONFLICT == 1
4793      if (numberConflictCuts > lastNumberConflictCuts) {
4794        double length = lengthConflictCuts / numberConflictCuts;
4795        printf("%d new conflict cuts - total %d - average length %g\n",
4796          numberConflictCuts - lastNumberConflictCuts,
4797          numberConflictCuts, length);
4798        lastNumberConflictCuts = numberConflictCuts;
4799      }
4800#endif
4801      if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4802        eventHappened_ = true; // exit
4803      }
4804    }
4805    // See if can stop on gap
4806    if (canStopOnGap()) {
4807      stoppedOnGap_ = true;
4808    }
4809
4810#ifdef CHECK_NODE_FULL
4811    verifyTreeNodes(tree_, *this);
4812#endif
4813#ifdef CHECK_CUT_COUNTS
4814    verifyCutCounts(tree_, *this);
4815#endif
4816    /*
4817          Now we come to the meat of the loop. To create the active subproblem, we'll
4818          pop the most promising node in the live set, rebuild the subproblem it
4819          represents, and then execute the current arm of the branch to create the
4820          active subproblem.
4821        */
4822    CbcNode *node = NULL;
4823#ifdef CBC_THREAD
4824    if (!parallelMode() || parallelMode() == -1) {
4825#endif
4826      node = tree_->bestNode(cutoff);
4827      // Possible one on tree worse than cutoff
4828      // Weird comparison function can leave ineligible nodes on tree
4829      if (!node || node->objectiveValue() > cutoff)
4830        continue;
4831      // Do main work of solving node here
4832      doOneNode(this, node, createdNode);
4833#ifdef JJF_ZERO
4834      if (node) {
4835        if (createdNode) {
4836          printf("Node %d depth %d, created %d depth %d\n",
4837            node->nodeNumber(), node->depth(),
4838            createdNode->nodeNumber(), createdNode->depth());
4839        } else {
4840          printf("Node %d depth %d,  no created node\n",
4841            node->nodeNumber(), node->depth());
4842        }
4843      } else if (createdNode) {
4844        printf("Node exhausted, created %d depth %d\n",
4845          createdNode->nodeNumber(), createdNode->depth());
4846      } else {
4847        printf("Node exhausted,  no created node\n");
4848        numberConsecutiveInfeasible = 2;
4849      }
4850#endif
4851      //if (createdNode)
4852      //numberConsecutiveInfeasible=0;
4853      //else
4854      //numberConsecutiveInfeasible++;
4855#ifdef CBC_THREAD
4856    } else if (parallelMode() > 0) {
4857      //lockThread();
4858      //node = tree_->bestNode(cutoff) ;
4859      // Possible one on tree worse than cutoff
4860      if (true || !node || node->objectiveValue() > cutoff) {
4861        assert(master_);
4862        if (master_) {
4863          int anyLeft = master_->waitForThreadsInTree(1);
4864          // may need to go round again
4865          if (anyLeft) {
4866            continue;
4867          } else {
4868            master_->stopThreads(-1);
4869          }
4870        }
4871      }
4872      //unlockThread();
4873    } else {
4874      // Deterministic parallel
4875      if ((tree_->size() < CoinMax(numberThreads_, 8) || hotstartSolution_) && !goneParallel) {
4876        node = tree_->bestNode(cutoff);
4877        // Possible one on tree worse than cutoff
4878        if (!node || node->objectiveValue() > cutoff)
4879          continue;
4880        // Do main work of solving node here
4881        doOneNode(this, node, createdNode);
4882        assert(createdNode);
4883        if (!createdNode->active()) {
4884          delete createdNode;
4885          createdNode = NULL;
4886        } else {
4887          // Say one more pointing to this
4888          node->nodeInfo()->increment();
4889          tree_->push(createdNode);
4890        }
4891        if (node->active()) {
4892          assert(node->nodeInfo());
4893          if (node->nodeInfo()->numberBranchesLeft()) {
4894            tree_->push(node);
4895          } else {
4896            node->setActive(false);
4897          }
4898        } else {
4899          if (node->nodeInfo()) {
4900            if (!node->nodeInfo()->numberBranchesLeft())
4901              node->nodeInfo()->allBranchesGone(); // can clean up
4902            // So will delete underlying stuff
4903            node->setActive(true);
4904          }
4905          delNode[nDeleteNode++] = node;
4906          node = NULL;
4907        }
4908        if (nDeleteNode >= MAX_DEL_NODE) {
4909          for (int i = 0; i < nDeleteNode; i++) {
4910            //printf("trying to del %d %x\n",i,delNode[i]);
4911            delete delNode[i];
4912            //printf("done to del %d %x\n",i,delNode[i]);
4913          }
4914          nDeleteNode = 0;
4915        }
4916      } else {
4917        // Split and solve
4918        master_->deterministicParallel();
4919        goneParallel = true;
4920      }
4921    }
4922#endif
4923  }
4924  if (nDeleteNode) {
4925    for (int i = 0; i < nDeleteNode; i++) {
4926      delete delNode[i];
4927    }
4928    nDeleteNode = 0;
4929  }
4930#ifdef CBC_THREAD
4931  if (master_) {
4932    master_->stopThreads(-1);
4933    master_->waitForThreadsInTree(2);
4934    // adjust time to allow for children on some systems
4935    //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4936  }
4937#endif
4938  /*
4939      End of the non-abort actions. The next block of code is executed if we've
4940      aborted because we hit one of the limits. Clean up by deleting the live set
4941      and break out of the node processing loop. Note that on an abort, node may
4942      have been pushed back onto the tree for further processing, in which case
4943      it'll be deleted in cleanTree. We need to check.
4944    */
4945  if (!(numberNodes_ < intParam_[CbcMaxNumNode] && numberSolutions_ < intParam_[CbcMaxNumSol] && !maximumSecondsReached() && !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))) {
4946    if (tree_->size()) {
4947      double dummyBest;
4948      tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4949#if 0 // Does not seem to be needed def CBC_THREAD
4950            if (parallelMode() > 0 && master_) {
4951              // see if any dangling nodes
4952              int numberThreads = master_->numberThreads();
4953              for (int i=0;i<numberThreads;i++) {
4954                CbcThread * child = master_->child(i);
4955                //if (child->createdNode())
4956                //printf("CHILD_NODE %p\n",child->createdNode());
4957                delete child->createdNode();
4958              }
4959            }
4960#endif
4961    }
4962    delete nextRowCut_;
4963    /* order is important here:
4964         * maximumSecondsReached() should be checked before eventHappened_ and
4965         * isNodeLimitReached() should be checked after eventHappened_
4966         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4967         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4968         */
4969    if (stoppedOnGap_) {
4970      messageHandler()->message(CBC_GAP, messages())
4971        << bestObjective_ - bestPossibleObjective_
4972        << dblParam_[CbcAllowableGap]
4973        << dblParam_[CbcAllowableFractionGap] * 100.0
4974        << CoinMessageEol;
4975      secondaryStatus_ = 2;
4976      status_ = 0;
4977    } else if (maximumSecondsReached()) {
4978      handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol;
4979      secondaryStatus_ = 4;
4980      status_ = 1;
4981    } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4982      handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol;
4983      secondaryStatus_ = 6;
4984      status_ = 1;
4985    } else if (isNodeLimitReached()) {
4986      handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol;
4987      secondaryStatus_ = 3;
4988      status_ = 1;
4989    } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4990      handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol;
4991      secondaryStatus_ = 8;
4992      status_ = 1;
4993    } else {
4994      handler_->message(CBC_EVENT, messages_) << CoinMessageEol;
4995      secondaryStatus_ = 5;
4996      status_ = 5;
4997    }
4998  }
4999#ifdef CBC_THREAD
5000  if (master_) {
5001    delete master_;
5002    master_ = NULL;
5003    masterThread_ = NULL;
5004  }
5005#endif
5006  /*
5007      That's it, we've exhausted the search tree, or broken out of the loop because
5008      we hit some limit on evaluation.
5009
5010      We may have got an intelligent tree so give it one more chance
5011    */
5012  // Tell solver we are not in Branch and Cut
5013  solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL);
5014  tree_->endSearch();
5015  //  If we did any sub trees - did we give up on any?
5016  if (numberStoppedSubTrees_)
5017    status_ = 1;
5018  numberNodes_ += numberExtraNodes_;
5019  numberIterations_ += numberExtraIterations_;
5020  if (eventHandler) {
5021    eventHandler->event(CbcEventHandler::endSearch);
5022  }
5023  if (!status_) {
5024    // Set best possible unless stopped on gap
5025    if (secondaryStatus_ != 2)
5026      bestPossibleObjective_ = bestObjective_;
5027    handler_->message(CBC_END_GOOD, messages_)
5028      << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
5029      << CoinMessageEol;
5030  } else {
5031    handler_->message(CBC_END, messages_)
5032      << bestObjective_ << bestPossibleObjective_
5033      << numberIterations_ << numberNodes_ << getCurrentSeconds()
5034      << CoinMessageEol;
5035  }
5036  if ((moreSpecialOptions_ & 4194304) != 0) {
5037    // Conflict cuts
5038    int numberCuts = globalCuts_.sizeRowCuts();
5039    int nConflict = 0;
5040    double sizeConflict = 0.0;
5041    for (int i = 0; i < numberCuts; i++) {
5042      OsiRowCut2 *cut = globalCuts_.cut(i);
5043      if (cut->whichRow() == 1) {
5044        nConflict++;
5045        sizeConflict += cut->row().getNumElements();
5046      }
5047    }
5048    if (nConflict) {
5049      sizeConflict /= nConflict;
5050      char general[200];
5051      sprintf(general, "%d conflict cuts generated - average length %g",
5052        nConflict, sizeConflict);
5053      messageHandler()->message(CBC_GENERAL,
5054        messages())
5055        << general << CoinMessageEol;
5056    }
5057  }
5058  if (numberStrongIterations_)
5059    handler_->message(CBC_STRONG_STATS, messages_)
5060      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
5061      << strongInfo_[1] << CoinMessageEol;
5062  if (!numberExtraNodes_)
5063    handler_->message(CBC_OTHER_STATS, messages_)
5064      << maximumDepthActual_
5065      << numberDJFixed_ << CoinMessageEol;
5066  else
5067    handler_->message(CBC_OTHER_STATS2, messages_)
5068      << maximumDepthActual_
5069      << numberDJFixed_ << numberFathoms_ << numberExtraNodes_ << numberExtraIterations_
5070      << CoinMessageEol;
5071#ifdef COIN_HAS_NTY
5072  if (symmetryInfo_)
5073    symmetryInfo_->statsOrbits(this, 1);
5074#endif
5075  if (doStatistics == 100) {
5076    for (int i = 0; i < numberObjects_; i++) {
5077      CbcSimpleIntegerDynamicPseudoCost *obj = dynamic_cast<CbcSimpleIntegerDynamicPseudoCost *>(object_[i]);
5078      if (obj)
5079        obj->print();
5080    }
5081  }
5082  if (statistics_) {
5083    // report in some way
5084    int *lookup = new int[numberObjects_];
5085    int i;
5086    for (i = 0; i < numberObjects_; i++)
5087      lookup[i] = -1;
5088    bool goodIds = false; //true;
5089    for (i = 0; i < numberObjects_; i++) {
5090      int iColumn = object_[i]->columnNumber();
5091      if (iColumn >= 0 && iColumn < numberColumns) {
5092        if (lookup[i] == -1) {
5093          lookup[i] = iColumn;
5094        } else {
5095          goodIds = false;
5096          break;
5097        }
5098      } else {
5099        goodIds = false;
5100        break;
5101      }
5102    }
5103    if (!goodIds) {
5104      delete[] lookup;
5105      lookup = NULL;
5106    }
5107    if (doStatistics >= 3) {
5108      printf("  node parent depth column   value                    obj      inf\n");
5109      for (i = 0; i < numberNodes2_; i++) {
5110        statistics_[i]->print(lookup);
5111      }
5112    }
5113    if (doStatistics > 1) {
5114      // Find last solution
5115      int k;
5116      for (k = numberNodes2_ - 1; k >= 0; k--) {
5117        if (statistics_[k]->endingObjective() != COIN_DBL_MAX && !statistics_[k]->endingInfeasibility())
5118          break;
5119      }
5120      if (k >= 0) {
5121        int depth = statistics_[k]->depth();
5122        int *which = new int[depth + 1];
5123        for (i = depth; i >= 0; i--) {
5124          which[i] = k;
5125          k = statistics_[k]->parentNode();
5126        }
5127        printf("  node parent depth column   value                    obj      inf\n");
5128        for (i = 0; i <= depth; i++) {
5129          statistics_[which[i]]->print(lookup);
5130        }
5131        delete[] which;
5132      }
5133    }
5134    // now summary
5135    int maxDepth = 0;
5136    double averageSolutionDepth = 0.0;
5137    int numberSolutions = 0;
5138    double averageCutoffDepth = 0.0;
5139    double averageSolvedDepth = 0.0;
5140    int numberCutoff = 0;
5141    int numberDown = 0;
5142    int numberFirstDown = 0;
5143    double averageInfDown = 0.0;
5144    double averageObjDown = 0.0;
5145    int numberCutoffDown = 0;
5146    int numberUp = 0;
5147    int numberFirstUp = 0;
5148    double averageInfUp = 0.0;
5149    double averageObjUp = 0.0;
5150    int numberCutoffUp = 0;
5151    double averageNumberIterations1 = 0.0;
5152    double averageValue = 0.0;
5153    for (i = 0; i < numberNodes2_; i++) {
5154      int depth = statistics_[i]->depth();
5155      int way = statistics_[i]->way();
5156      double value = statistics_[i]->value();
5157      double startingObjective = statistics_[i]->startingObjective();
5158      int startingInfeasibility = statistics_[i]->startingInfeasibility();
5159      double endingObjective = statistics_[i]->endingObjective();
5160      int endingInfeasibility = statistics_[i]->endingInfeasibility();
5161      maxDepth = CoinMax(depth, maxDepth);
5162      // Only for completed
5163      averageNumberIterations1 += statistics_[i]->numberIterations();
5164      averageValue += value;
5165      if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
5166        numberSolutions++;
5167        averageSolutionDepth += depth;
5168      }
5169      if (endingObjective == COIN_DBL_MAX) {
5170        numberCutoff++;
5171        averageCutoffDepth += depth;
5172        if (way < 0) {
5173          numberDown++;
5174          numberCutoffDown++;
5175          if (way == -1)
5176            numberFirstDown++;
5177        } else {
5178          numberUp++;
5179          numberCutoffUp++;
5180          if (way == 1)
5181            numberFirstUp++;
5182        }
5183      } else {
5184        averageSolvedDepth += depth;
5185        if (way < 0) {
5186          numberDown++;
5187          averageInfDown += startingInfeasibility - endingInfeasibility;
5188          averageObjDown += endingObjective - startingObjective;
5189          if (way == -1)
5190            numberFirstDown++;
5191        } else {
5192          numberUp++;
5193          averageInfUp += startingInfeasibility - endingInfeasibility;
5194          averageObjUp += endingObjective - startingObjective;
5195          if (way == 1)
5196            numberFirstUp++;
5197        }
5198      }
5199    }
5200    // Now print
5201    if (numberSolutions)
5202      averageSolutionDepth /= static_cast<double>(numberSolutions);
5203    int numberSolved = numberNodes2_ - numberCutoff;
5204    double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
5205      - numberIterationsAtContinuous;
5206    if (numberCutoff) {
5207      averageCutoffDepth /= static_cast<double>(numberCutoff);
5208      averageNumberIterations2 /= static_cast<double>(numberCutoff);
5209    }
5210    if (numberNodes2_)
5211      averageValue /= static_cast<double>(numberNodes2_);
5212    if (numberSolved) {
5213      averageNumberIterations1 /= static_cast<double>(numberSolved);
5214      averageSolvedDepth /= static_cast<double>(numberSolved);
5215    }
5216    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
5217      numberSolutions, averageSolutionDepth);
5218    printf("average value of variable being branched on was %g\n",
5219      averageValue);
5220    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
5221      numberCutoff, averageCutoffDepth, averageNumberIterations2);
5222    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
5223      numberSolved, averageSolvedDepth, averageNumberIterations1);
5224    if (numberDown) {
5225      averageInfDown /= static_cast<double>(numberDown);
5226      averageObjDown /= static_cast<double>(numberDown);
5227    }
5228    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5229      numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
5230      averageInfDown, averageObjDown);
5231    if (numberUp) {
5232      averageInfUp /= static_cast<double>(numberUp);
5233      averageObjUp /= static_cast<double>(numberUp);
5234    }
5235    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5236      numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
5237      averageInfUp, averageObjUp);
5238    for (i = 0; i < numberNodes2_; i++)
5239      delete statistics_[i];
5240    delete[] statistics_;
5241    statistics_ = NULL;
5242    maximumStatistics_ = 0;
5243    delete[] lookup;
5244  }
5245  /*
5246      If we think we have a solution, restore and confirm it with a call to
5247      setBestSolution().  We need to reset the cutoff value so as not to fathom
5248      the solution on bounds.  Note that calling setBestSolution( ..., true)
5249      leaves the continuousSolver_ bounds vectors fixed at the solution value.
5250
5251      Running resolve() here is a failsafe --- setBestSolution has already
5252      reoptimised using the continuousSolver_. If for some reason we fail to
5253      prove optimality, run the problem again after instructing the solver to
5254      tell us more.
5255
5256      If all looks good, replace solver_ with continuousSolver_, so that the
5257      outside world will be able to obtain information about the solution using
5258      public methods.
5259
5260      Don't replace if we are trying to save cuts
5261    */
5262  if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && ((specialOptions_ & 8388608) == 0 || (specialOptions_ & 2048) != 0)) {
5263    setCutoff(1.0e50); // As best solution should be worse than cutoff
5264    // change cutoff as constraint if wanted
5265    if (cutoffRowNumber_ >= 0) {
5266      if (solver_->getNumRows() > cutoffRowNumber_)
5267        solver_->setRowUpper(cutoffRowNumber_, 1.0e50);
5268    }
5269    // also in continuousSolver_
5270    if (continuousSolver_) {
5271      // Solvers know about direction
5272      double direction = solver_->getObjSense();
5273      continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
5274    }
5275    phase_ = 5;
5276    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
5277    if ((specialOptions_ & 4) == 0)
5278      bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
5279    setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1);
5280    continuousSolver_->resolve();
5281    // Deal with funny variables
5282    if ((moreSpecialOptions2_ & 32768) != 0)
5283      cleanBounds(continuousSolver_, NULL);
5284    if (!continuousSolver_->isProvenOptimal()) {
5285      continuousSolver_->messageHandler()->setLogLevel(2);
5286      continuousSolver_->initialSolve();
5287    }
5288    delete solver_;
5289    // above deletes solverCharacteristics_
5290    solverCharacteristics_ = NULL;
5291    solver_ = continuousSolver_;
5292    setPointers(solver_);
5293    continuousSolver_ = NULL;
5294  }
5295  /*
5296      Clean up dangling objects. continuousSolver_ may already be toast.
5297    */
5298  delete lastws;
5299  if (saveObjects) {
5300    for (int i = 0; i < numberObjects_; i++)
5301      delete saveObjects[i];
5302    delete[] saveObjects;
5303  }
5304  numberStrong_ = saveNumberStrong;
5305  numberBeforeTrust_ = saveNumberBeforeTrust;
5306  delete[] whichGenerator_;
5307  whichGenerator_ = NULL;
5308  delete[] lowerBefore;
5309  delete[] upperBefore;
5310  delete[] walkback_;
5311  walkback_ = NULL;
5312  delete[] lastNodeInfo_;
5313  lastNodeInfo_ = NULL;
5314  delete[] lastNumberCuts_;
5315  lastNumberCuts_ = NULL;
5316  delete[] lastCut_;
5317  lastCut_ = NULL;
5318  delete[] addedCuts_;
5319  addedCuts_ = NULL;
5320  //delete persistentInfo;
5321  // Get rid of characteristics
5322  solverCharacteristics_ = NULL;
5323  if (continuousSolver_) {
5324    delete continuousSolver_;
5325    continuousSolver_ = NULL;
5326  }
5327  /*
5328      Destroy global cuts by replacing with an empty OsiCuts object.
5329    */
5330  globalCuts_ = CbcRowCuts();
5331  delete globalConflictCuts_;
5332  globalConflictCuts_ = NULL;
5333  if (!bestSolution_ && (specialOptions_ & 8388608) == 0 && false) {
5334    // make sure lp solver is infeasible
5335    int numberColumns = solver_->getNumCols();
5336    const double *columnLower = solver_->getColLower();
5337    int iColumn;
5338    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5339      if (solver_->isInteger(iColumn))
5340        solver_->setColUpper(iColumn, columnLower[iColumn]);
5341    }
5342    solver_->initialSolve();
5343  }
5344#ifdef COIN_HAS_CLP
5345  {
5346    OsiClpSolverInterface *clpSolver
5347      = dynamic_cast<OsiClpSolverInterface *>(solver_);
5348    if (clpSolver) {
5349      // Possible restore of pivot method
5350      if (savePivotMethod) {
5351        // model may have changed
5352        savePivotMethod->setModel(NULL);
5353        clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
5354        delete savePivotMethod;
5355      }
5356      clpSolver->setLargestAway(-1.0);
5357    }
5358  }
5359#endif
5360  if ((fastNodeDepth_ >= 1000000 || (moreSpecialOptions_ & 33554432) != 0)
5361    && !parentModel_) {
5362    // delete object off end
5363    delete object_[numberObjects_];
5364    if ((moreSpecialOptions_ & 33554432) == 0)
5365      fastNodeDepth_ -= 1000000;
5366  }
5367  delete saveSolver;
5368  // Undo preprocessing performed during BaB.
5369  if (strategy_ && strategy_->preProcessState() > 0) {
5370    // undo preprocessing
5371    CglPreProcess *process = strategy_->process();
5372    assert(process);
5373    int n = originalSolver->getNumCols();
5374    if (bestSolution_) {
5375      delete[] bestSolution_;
5376      bestSolution_ = new double[n];
5377      process->postProcess(*solver_);
5378    }
5379    strategy_->deletePreProcess();
5380    // Solution now back in originalSolver
5381    delete solver_;
5382    solver_ = originalSolver;
5383    if (bestSolution_) {
5384      bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
5385      memcpy(bestSolution_, solver_->getColSolution(), n * sizeof(double));
5386    }
5387    // put back original objects if there were any
5388    if (originalObject) {
5389      int iColumn;
5390      assert(ownObjects_);
5391      for (iColumn = 0; iColumn < numberObjects_; iColumn++)
5392        delete object_[iColumn];
5393      delete[] object_;
5394      numberObjects_ = numberOriginalObjects;
5395      object_ = originalObject;
5396      delete[] integerVariable_;
5397      numberIntegers_ = 0;
5398      for (iColumn = 0; iColumn < n; iColumn++) {
5399        if (solver_->isInteger(iColumn))
5400          numberIntegers_++;
5401      }
5402      integerVariable_ = new int[numberIntegers_];
5403      numberIntegers_ = 0;
5404      for (iColumn = 0; iColumn < n; iColumn++) {
5405        if (solver_->isInteger(iColumn))
5406          integerVariable_[numberIntegers_++] = iColumn;
5407      }
5408    }
5409  }
5410  if (flipObjective)
5411    flipModel();
5412#ifdef COIN_HAS_CLP
5413  {
5414    OsiClpSolverInterface *clpSolver
5415      = dynamic_cast<OsiClpSolverInterface *>(solver_);
5416    if (clpSolver)
5417      clpSolver->setFakeObjective(reinterpret_cast<double *>(NULL));
5418  }
5419#endif
5420  moreSpecialOptions_ = saveMoreSpecialOptions;
5421  return;
5422}
5423
5424// Solve the initial LP relaxation
5425void CbcModel::initialSolve()
5426{
5427  assert(solver_);
5428  // Double check optimization directions line up
5429  dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
5430  // Check if bounds are all integral (as may get messed up later)
5431  checkModel();
5432  if (!solverCharacteristics_) {
5433    OsiBabSolver *solverCharacteristics = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
5434    if (solverCharacteristics) {
5435      solverCharacteristics_ = solverCharacteristics;
5436    } else {
5437      // replace in solver
5438      OsiBabSolver defaultC;
5439      solver_->setAuxiliaryInfo(&defaultC);
5440      solverCharacteristics_ = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
5441    }
5442  }
5443  solverCharacteristics_->setSolver(solver_);
5444  solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL);
5445  solver_->initialSolve();
5446  solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL);
5447  if (!solver_->isProvenOptimal())
5448    solver_->resolve();
5449  // But set up so Jon Lee will be happy
5450  status_ = -1;
5451  secondaryStatus_ = -1;
5452  originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
5453  bestPossibleObjective_ = originalContinuousObjective_;
5454  if (solver_->isProvenDualInfeasible())
5455    originalContinuousObjective_ = -COIN_DBL_MAX;