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

Last change on this file since 2491 was 2491, checked in by forrest, 21 months ago

better SOS in mipstart, ctrl-c back, improve symmetric

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 703.3 KB
Line 
1/* $Id: CbcModel.cpp 2491 2019-02-12 13:05:24Z forrest $ */
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      if ( this->keepNamesPreproc == false )
2000        solver_->setIntParam( OsiNameDiscipline, 0 );
2001      ClpSimplex *clpSimplex = clpSolver->getModelPtr();
2002      if ((specialOptions_ & 32) == 0) {
2003        // take off names (unless going to be saving)
2004        int nameDisc; solver_->getIntParam( OsiNameDiscipline, nameDisc );
2005        if ( (numberAnalyzeIterations_ >= 0 || (-numberAnalyzeIterations_ & 64) == 0) && (!nameDisc) )
2006          clpSimplex->dropNames();
2007      }
2008      // no crunch if mostly continuous
2009      if ((clpSolver->specialOptions() & (1 + 8)) != (1 + 8)) {
2010        int numberColumns = solver_->getNumCols();
2011        if (numberColumns > 1000 && numberIntegers_ * 4 < numberColumns)
2012          clpSolver->setSpecialOptions(clpSolver->specialOptions() & (~1));
2013      }
2014      //#define NO_CRUNCH
2015#ifdef NO_CRUNCH
2016      printf("TEMP switching off crunch\n");
2017      int iOpt = clpSolver->specialOptions();
2018      iOpt &= ~1;
2019      iOpt |= 65536;
2020      clpSolver->setSpecialOptions(iOpt);
2021#endif
2022    }
2023  }
2024#endif
2025  bool feasible;
2026  numberSolves_ = 0;
2027  {
2028    // check
2029    int numberOdd = 0;
2030    for (int i = 0; i < numberObjects_; i++) {
2031      CbcSimpleInteger *obj = dynamic_cast< CbcSimpleInteger * >(object_[i]);
2032      if (!obj)
2033        numberOdd++;
2034    }
2035    if (numberOdd)
2036      moreSpecialOptions_ |= 1073741824;
2037  }
2038  // If NLP then we assume already solved outside branchAndbound
2039  if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
2040    feasible = resolve(NULL, 0) != 0;
2041  } else {
2042    // pick up given status
2043    feasible = (solver_->isProvenOptimal() && !solver_->isDualObjectiveLimitReached());
2044  }
2045  if (problemFeasibility_->feasible(this, 0) < 0) {
2046    feasible = false; // pretend infeasible
2047  }
2048  numberSavedSolutions_ = 0;
2049  int saveNumberStrong = numberStrong_;
2050  int saveNumberBeforeTrust = numberBeforeTrust_;
2051  /*
2052      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2053      continue with processing the root node.
2054    */
2055  if (!feasible) {
2056    status_ = 0;
2057    if (!solver_->isProvenDualInfeasible()) {
2058      handler_->message(CBC_INFEAS, messages_) << CoinMessageEol;
2059      secondaryStatus_ = 1;
2060    } else {
2061      handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol;
2062      secondaryStatus_ = 7;
2063    }
2064    originalContinuousObjective_ = COIN_DBL_MAX;
2065    if (bestSolution_ && ((specialOptions_ & 8388608) == 0 || (specialOptions_ & 2048) != 0)) {
2066      // best solution found by various heuristics - set solution
2067      char general[200];
2068      sprintf(general, "Solution of %g already found by heuristic",
2069        bestObjective_);
2070      messageHandler()->message(CBC_GENERAL,
2071        messages())
2072        << general << CoinMessageEol;
2073      setCutoff(1.0e50); // As best solution should be worse than cutoff
2074      // change cutoff as constraint if wanted
2075      if (cutoffRowNumber_ >= 0) {
2076        if (solver_->getNumRows() > cutoffRowNumber_)
2077          solver_->setRowUpper(cutoffRowNumber_, 1.0e50);
2078      }
2079      // also in continuousSolver_
2080      if (continuousSolver_) {
2081        // Solvers know about direction
2082        double direction = solver_->getObjSense();
2083        continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
2084      } else {
2085        continuousSolver_ = solver_->clone();
2086      }
2087      phase_ = 5;
2088      double increment = getDblParam(CbcModel::CbcCutoffIncrement);
2089      if ((specialOptions_ & 4) == 0)
2090        bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
2091      setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1);
2092      continuousSolver_->resolve();
2093      if (!continuousSolver_->isProvenOptimal()) {
2094        continuousSolver_->messageHandler()->setLogLevel(2);
2095        continuousSolver_->initialSolve();
2096      }
2097      delete solver_;
2098      solverCharacteristics_ = NULL;
2099      solver_ = continuousSolver_;
2100      setPointers(solver_);
2101      continuousSolver_ = NULL;
2102    }
2103    solverCharacteristics_ = NULL;
2104    if (flipObjective)
2105      flipModel();
2106    return;
2107  } else if (!numberObjects_) {
2108    // nothing to do
2109    // Undo preprocessing performed during BaB.
2110    if (strategy_ && strategy_->preProcessState() > 0) {
2111      // undo preprocessing
2112      CglPreProcess *process = strategy_->process();
2113      assert(process);
2114      int n = originalSolver->getNumCols();
2115      if (bestSolution_) {
2116        delete[] bestSolution_;
2117        bestSolution_ = new double[n];
2118        process->postProcess(*solver_);
2119      }
2120      strategy_->deletePreProcess();
2121      // Solution now back in originalSolver
2122      delete solver_;
2123      solver_ = originalSolver;
2124      if (bestSolution_) {
2125        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2126        memcpy(bestSolution_, solver_->getColSolution(), n * sizeof(double));
2127      }
2128      // put back original objects if there were any
2129      if (originalObject) {
2130        int iColumn;
2131        assert(ownObjects_);
2132        for (iColumn = 0; iColumn < numberObjects_; iColumn++)
2133          delete object_[iColumn];
2134        delete[] object_;
2135        numberObjects_ = numberOriginalObjects;
2136        object_ = originalObject;
2137        delete[] integerVariable_;
2138        numberIntegers_ = 0;
2139        for (iColumn = 0; iColumn < n; iColumn++) {
2140          if (solver_->isInteger(iColumn))
2141            numberIntegers_++;
2142        }
2143        integerVariable_ = new int[numberIntegers_];
2144        numberIntegers_ = 0;
2145        for (iColumn = 0; iColumn < n; iColumn++) {
2146          if (solver_->isInteger(iColumn))
2147            integerVariable_[numberIntegers_++] = iColumn;
2148        }
2149      }
2150    }
2151    if (flipObjective)
2152      flipModel();
2153    solverCharacteristics_ = NULL;
2154    bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2155    int numberColumns = solver_->getNumCols();
2156    delete[] bestSolution_;
2157    bestSolution_ = new double[numberColumns];
2158    CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2159    return;
2160  }
2161  /*
2162      See if we're using the Osi side of the branching hierarchy. If so, either
2163      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2164      first case, CbcModel owns the objects on the object_ list. In the second
2165      case, the solver holds the objects and object_ simply points to the
2166      solver's list.
2167
2168      080417 The conversion code here (the block protected by `if (obj)') cannot
2169      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2170      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2171      -> all other Cbc object classes. It's structurally impossible for any Osi
2172      object to descend from CbcObject. The only thing I can see is that this is
2173      really dead code, and object detection is now handled from the Osi side.
2174    */
2175  // Convert to Osi if wanted
2176  //OsiBranchingInformation * persistentInfo = NULL;
2177  if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2178    //persistentInfo = new OsiBranchingInformation(solver_);
2179    if (numberOriginalObjects) {
2180      for (int iObject = 0; iObject < numberObjects_; iObject++) {
2181        CbcObject *obj = dynamic_cast< CbcObject * >(object_[iObject]);
2182        if (obj) {
2183          CbcSimpleInteger *obj2 = dynamic_cast< CbcSimpleInteger * >(obj);
2184          if (obj2) {
2185            // back to Osi land
2186            object_[iObject] = obj2->osiObject();
2187            delete obj;
2188          } else {
2189            OsiSimpleInteger *obj3 = dynamic_cast< OsiSimpleInteger * >(obj);
2190            if (!obj3) {
2191              OsiSOS *obj4 = dynamic_cast< OsiSOS * >(obj);
2192              if (!obj4) {
2193                CbcSOS *obj5 = dynamic_cast< CbcSOS * >(obj);
2194                if (obj5) {
2195                  // back to Osi land
2196                  object_[iObject] = obj5->osiObject(solver_);
2197                } else {
2198                  printf("Code up CbcObject type in Osi land\n");
2199                  abort();
2200                }
2201              }
2202            }
2203          }
2204        }
2205      }
2206      // and add to solver
2207      //if (!solver_->numberObjects()) {
2208      solver_->addObjects(numberObjects_, object_);
2209      //} else {
2210      //if (solver_->numberObjects()!=numberOriginalObjects) {
2211      //printf("should have trapped that solver has objects before\n");
2212      //abort();
2213      //}
2214      //}
2215    } else {
2216      /*
2217              As of 080104, findIntegersAndSOS is misleading --- the default OSI
2218              implementation finds only integers.
2219            */
2220      // do from solver
2221      deleteObjects(false);
2222      solver_->findIntegersAndSOS(false);
2223      numberObjects_ = solver_->numberObjects();
2224      object_ = solver_->objects();
2225      ownObjects_ = false;
2226    }
2227    branchingMethod_->chooseMethod()->setSolver(solver_);
2228  }
2229  // take off heuristics if have to (some do not work with SOS, for example)
2230  // object should know what's safe.
2231  {
2232    int numberOdd = 0;
2233    int numberSOS = 0;
2234    for (int i = 0; i < numberObjects_; i++) {
2235      if (!object_[i]->canDoHeuristics())
2236        numberOdd++;
2237      CbcSOS *obj = dynamic_cast< CbcSOS * >(object_[i]);
2238      if (obj)
2239        numberSOS++;
2240    }
2241    if (numberOdd) {
2242      if (numberHeuristics_ && (specialOptions_ & 1024) == 0) {
2243        int k = 0;
2244        for (int i = 0; i < numberHeuristics_; i++) {
2245          if (!heuristic_[i]->canDealWithOdd())
2246            delete heuristic_[i];
2247          else
2248            heuristic_[k++] = heuristic_[i];
2249        }
2250        if (!k) {
2251          delete[] heuristic_;
2252          heuristic_ = NULL;
2253        }
2254        numberHeuristics_ = k;
2255        handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol;
2256      }
2257      // If odd switch off AddIntegers
2258      specialOptions_ &= ~65536;
2259      // switch off fast nodes for now
2260      fastNodeDepth_ = -1;
2261      moreSpecialOptions_ &= ~33554432; // no diving
2262    } else if (numberSOS) {
2263      specialOptions_ |= 128; // say can do SOS in dynamic mode
2264      // switch off fast nodes for now
2265      fastNodeDepth_ = -1;
2266      moreSpecialOptions_ &= ~33554432; // no diving
2267    }
2268    if (numberThreads_ > 0) {
2269      /* switch off fast nodes for now
2270               Trouble is that by time mini bab finishes code is
2271               looking at a different node
2272             */
2273      fastNodeDepth_ = -1;
2274    }
2275  }
2276  // Save objective (just so user can access it)
2277  originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
2278  bestPossibleObjective_ = originalContinuousObjective_;
2279  sumChangeObjective1_ = 0.0;
2280  sumChangeObjective2_ = 0.0;
2281  /*
2282      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2283      Assuming it recognises the problem, when called upon it will check a cut to
2284      see if it cuts off the optimal answer.
2285    */
2286  // If debugger exists set specialOptions_ bit
2287  if (solver_->getRowCutDebuggerAlways()) {
2288    specialOptions_ |= 1;
2289  }
2290
2291#ifdef CBC_DEBUG
2292  if ((specialOptions_ & 1) == 0)
2293    solver_->activateRowCutDebugger(problemName.c_str());
2294  if (solver_->getRowCutDebuggerAlways())
2295    specialOptions_ |= 1;
2296#endif
2297
2298  /*
2299      Begin setup to process a feasible root node.
2300    */
2301  bestObjective_ = CoinMin(bestObjective_, 1.0e50);
2302  if (!bestSolution_) {
2303    numberSolutions_ = 0;
2304    numberHeuristicSolutions_ = 0;
2305  }
2306  stateOfSearch_ = 0;
2307  // Everything is minimization
2308  {
2309    // needed to sync cutoffs
2310    double value;
2311    solver_->getDblParam(OsiDualObjectiveLimit, value);
2312    dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2313  }
2314  double cutoff = getCutoff();
2315  double direction = solver_->getObjSense();
2316  dblParam_[CbcOptimizationDirection] = direction;
2317  if (cutoff < 1.0e20 && direction < 0.0)
2318    messageHandler()->message(CBC_CUTOFF_WARNING1,
2319      messages())
2320      << cutoff << -cutoff << CoinMessageEol;
2321  if (cutoff > bestObjective_)
2322    cutoff = bestObjective_;
2323  setCutoff(cutoff);
2324  /*
2325      We probably already have a current solution, but just in case ...
2326    */
2327  int numberColumns = getNumCols();
2328  if (!currentSolution_)
2329    currentSolution_ = new double[numberColumns];
2330  testSolution_ = currentSolution_;
2331  /*
2332      Create a copy of the solver, thus capturing the original (root node)
2333      constraint system (aka the continuous system).
2334    */
2335  delete continuousSolver_;
2336  continuousSolver_ = solver_->clone();
2337#ifdef CONFLICT_CUTS
2338  if ((moreSpecialOptions_ & 4194304) != 0) {
2339#ifdef COIN_HAS_CLP
2340    OsiClpSolverInterface *clpSolver
2341      = dynamic_cast< OsiClpSolverInterface * >(solver_);
2342    if (clpSolver) {
2343      int specialOptions = clpSolver->getModelPtr()->specialOptions();
2344      // 2097152 switches on rays in crunch
2345      if (!parentModel_)
2346        clpSolver->getModelPtr()->setSpecialOptions(specialOptions | 32 | 2097152);
2347      else
2348        clpSolver->getModelPtr()->setSpecialOptions(specialOptions & ~(32 | 2097152));
2349    }
2350  }
2351#endif
2352#endif
2353#ifdef COIN_HAS_NTY
2354  // maybe allow on fix and restart later
2355  if ((moreSpecialOptions2_ & (128 | 256)) != 0 && !parentModel_) {
2356    symmetryInfo_ = new CbcSymmetry();
2357    symmetryInfo_->setupSymmetry(this);
2358    int numberGenerators = symmetryInfo_->statsOrbits(this, 0);
2359    if (!symmetryInfo_->numberUsefulOrbits() && (moreSpecialOptions2_ & (128 | 256)) != (128 | 256)) {
2360      delete symmetryInfo_;
2361      symmetryInfo_ = NULL;
2362      moreSpecialOptions2_ &= ~(128 | 256);
2363    }
2364    if ((moreSpecialOptions2_ & (128 | 256)) == (128 | 256)) {
2365      //moreSpecialOptions2_ &= ~256;
2366    }
2367  }
2368#endif
2369
2370  // add cutoff as constraint if wanted
2371  if (cutoffRowNumber_ == -2) {
2372    if (!parentModel_) {
2373      int numberColumns = solver_->getNumCols();
2374      double *obj = CoinCopyOfArray(solver_->getObjCoefficients(), numberColumns);
2375      int *indices = new int[numberColumns];
2376      int n = 0;
2377      for (int i = 0; i < numberColumns; i++) {
2378        if (obj[i]) {
2379          indices[n] = i;
2380          obj[n++] = obj[i];
2381        }
2382      }
2383      if (n) {
2384        double cutoff = getCutoff();
2385        // relax a little bit
2386        cutoff += 1.0e-4;
2387        double offset;
2388        solver_->getDblParam(OsiObjOffset, offset);
2389        cutoffRowNumber_ = solver_->getNumRows();
2390        solver_->addRow(n, indices, obj, -COIN_DBL_MAX, CoinMin(cutoff, 1.0e25) + offset);
2391      } else {
2392        // no objective!
2393        cutoffRowNumber_ = -1;
2394      }
2395      delete[] indices;
2396      delete[] obj;
2397    } else {
2398      // switch off
2399      cutoffRowNumber_ = -1;
2400    }
2401  }
2402  numberRowsAtContinuous_ = getNumRows();
2403  solver_->saveBaseModel();
2404  /*
2405      Check the objective to see if we can deduce a nontrivial increment. If
2406      it's better than the current value for CbcCutoffIncrement, it'll be
2407      installed.
2408    */
2409  if (solverCharacteristics_->reducedCostsAccurate())
2410    analyzeObjective();
2411  {
2412    // may be able to change cutoff now
2413    double cutoff = getCutoff();
2414    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
2415    if (cutoff > bestObjective_ - increment) {
2416      cutoff = bestObjective_ - increment;
2417      setCutoff(cutoff);
2418    }
2419  }
2420#ifdef COIN_HAS_CLP
2421  // Possible save of pivot method
2422  ClpDualRowPivot *savePivotMethod = NULL;
2423  {
2424    // pass tolerance and increment to solver
2425    OsiClpSolverInterface *clpSolver
2426      = dynamic_cast< OsiClpSolverInterface * >(solver_);
2427    if (clpSolver)
2428      clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2429#ifdef CLP_RESOLVE
2430    if ((moreSpecialOptions_ & 1048576) != 0 && !parentModel_ && clpSolver) {
2431      resolveClp(clpSolver, 0);
2432    }
2433#endif
2434  }
2435#endif
2436  /*
2437      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2438      the active subproblem. whichGenerator will be used to record the generator
2439      that produced a given cut.
2440    */
2441#define INITIAL_MAXIMUM_WHICH 1000
2442  maximumWhich_ = INITIAL_MAXIMUM_WHICH;
2443  delete[] whichGenerator_;
2444  whichGenerator_ = new int[maximumWhich_];
2445  memset(whichGenerator_, 0, maximumWhich_ * sizeof(int));
2446  maximumNumberCuts_ = 0;
2447  currentNumberCuts_ = 0;
2448  delete[] addedCuts_;
2449  addedCuts_ = NULL;
2450  OsiObject **saveObjects = NULL;
2451  maximumRows_ = numberRowsAtContinuous_;
2452  currentDepth_ = 0;
2453  workingBasis_.resize(maximumRows_, numberColumns);
2454  /*
2455      Set up an empty heap and associated data structures to hold the live set
2456      (problems which require further exploration).
2457    */
2458  CbcCompareDefault *compareActual
2459    = dynamic_cast< CbcCompareDefault * >(nodeCompare_);
2460  if (compareActual) {
2461    compareActual->setBestPossible(direction * solver_->getObjValue());
2462    compareActual->setCutoff(getCutoff());
2463#ifdef JJF_ZERO
2464    if (false && !numberThreads_ && !parentModel_) {
2465      printf("CbcTreeArray ? threads ? parentArray\n");
2466      // Setup new style tree
2467      delete tree_;
2468      tree_ = new CbcTreeArray();
2469    }
2470#endif
2471  }
2472  tree_->setComparison(*nodeCompare_);
2473  /*
2474      Used to record the path from a node to the root of the search tree, so that
2475      we can then traverse from the root to the node when restoring a subproblem.
2476    */
2477  maximumDepth_ = 10;
2478  delete[] walkback_;
2479  walkback_ = new CbcNodeInfo *[maximumDepth_];
2480  lastDepth_ = 0;
2481  delete[] lastNodeInfo_;
2482  lastNodeInfo_ = new CbcNodeInfo *[maximumDepth_];
2483  delete[] lastNumberCuts_;
2484  lastNumberCuts_ = new int[maximumDepth_];
2485  maximumCuts_ = 100;
2486  lastNumberCuts2_ = 0;
2487  delete[] lastCut_;
2488  lastCut_ = new const OsiRowCut *[maximumCuts_];
2489  /*
2490      Used to generate bound edits for CbcPartialNodeInfo.
2491    */
2492  double *lowerBefore = new double[numberColumns];
2493  double *upperBefore = new double[numberColumns];
2494  /*
2495    Set up to run heuristics and generate cuts at the root node. The heavy
2496    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
2497
2498    To start, tell cut generators they can be a bit more aggressive at the
2499    root node.
2500
2501    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2502        with cuts at root'. Is phase_ = 1 the correct indication when
2503        doHeurisiticsAtRoot is called to run heuristics outside of the main
2504        cut / heurisitc / reoptimise loop in solveWithCuts?
2505
2506      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2507      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2508      fixing) until no cuts are generated, the change in objective falls off,  or
2509      the limit on the number of rounds of cut generation is exceeded.
2510
2511      At the end of all this, any cuts will be recorded in cuts and also
2512      installed in the solver's constraint system. We'll have reoptimised, and
2513      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2514      adjusted accordingly).
2515
2516      Tell cut generators they can be a bit more aggressive at root node
2517
2518      TODO: Why don't we make a copy of the solution after solveWithCuts?
2519      TODO: If numberUnsatisfied == 0, don't we have a solution?
2520    */
2521  phase_ = 1;
2522  int iCutGenerator;
2523  for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2524    // If parallel switch off global cuts
2525    if (numberThreads_) {
2526      generator_[iCutGenerator]->setGlobalCuts(false);
2527      generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2528    }
2529    CglCutGenerator *generator = generator_[iCutGenerator]->generator();
2530    generator->setAggressiveness(generator->getAggressiveness() + 100);
2531    if (!generator->canDoGlobalCuts())
2532      generator->setGlobalCuts(false);
2533  }
2534  OsiCuts cuts;
2535  int anyAction = -1;
2536  numberOldActiveCuts_ = 0;
2537  numberNewCuts_ = 0;
2538  // Array to mark solution
2539  delete[] usedInSolution_;
2540  usedInSolution_ = new int[numberColumns];
2541  CoinZeroN(usedInSolution_, numberColumns);
2542  /*
2543      For printing totals and for CbcNode (numberNodes_)
2544    */
2545  numberIterations_ = 0;
2546  numberNodes_ = 0;
2547  numberNodes2_ = 0;
2548  maximumStatistics_ = 0;
2549  maximumDepthActual_ = 0;
2550  numberDJFixed_ = 0.0;
2551  if (!parentModel_) {
2552    if ((specialOptions_ & 262144) != 0) {
2553      // create empty stored cuts
2554      //storedRowCuts_ = new CglStored(solver_->getNumCols());
2555    } else if ((specialOptions_ & 524288) != 0 && storedRowCuts_) {
2556      // tighten and set best solution
2557      // A) tight bounds on integer variables
2558      /*
2559                storedRowCuts_ are coming in from outside, probably for nonlinear.
2560              John was unsure about origin.
2561            */
2562      const double *lower = solver_->getColLower();
2563      const double *upper = solver_->getColUpper();
2564      const double *tightLower = storedRowCuts_->tightLower();
2565      const double *tightUpper = storedRowCuts_->tightUpper();
2566      int nTightened = 0;
2567      for (int i = 0; i < numberIntegers_; i++) {
2568        int iColumn = integerVariable_[i];
2569        if (tightLower[iColumn] > lower[iColumn]) {
2570          nTightened++;
2571          solver_->setColLower(iColumn, tightLower[iColumn]);
2572        }
2573        if (tightUpper[iColumn] < upper[iColumn]) {
2574          nTightened++;
2575          solver_->setColUpper(iColumn, tightUpper[iColumn]);
2576        }
2577      }
2578      if (nTightened)
2579        COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2580      if (storedRowCuts_->bestObjective() < bestObjective_) {
2581        // B) best solution
2582        double objValue = storedRowCuts_->bestObjective();
2583        setBestSolution(CBC_SOLUTION, objValue,
2584          storedRowCuts_->bestSolution());
2585        // Do heuristics
2586        // Allow RINS
2587        for (int i = 0; i < numberHeuristics_; i++) {
2588          CbcHeuristicRINS *rins
2589            = dynamic_cast< CbcHeuristicRINS * >(heuristic_[i]);
2590          if (rins) {
2591            rins->setLastNode(-100);
2592          }
2593        }
2594      }
2595    }
2596  }
2597#ifdef SWITCH_VARIABLES
2598  // see if any switching variables
2599  if (numberIntegers_ < solver_->getNumCols())
2600    findSwitching();
2601#endif
2602  /*
2603      Run heuristics at the root. This is the only opportunity to run FPump; it
2604      will be removed from the heuristics list by doHeuristicsAtRoot.
2605    */
2606  // See if multiple runs wanted
2607  CbcModel **rootModels = NULL;
2608  if (!parentModel_ && multipleRootTries_ % 100) {
2609    double rootTimeCpu = CoinCpuTime();
2610    double startTimeRoot = CoinGetTimeOfDay();
2611    int numberRootThreads = 1;
2612    /* undocumented fine tuning
2613         aabbcc where cc is number of tries
2614         bb if nonzero is number of threads
2615         aa if nonzero just do heuristics
2616      */
2617    int numberModels = multipleRootTries_ % 100;
2618#ifdef CBC_THREAD
2619    numberRootThreads = (multipleRootTries_ / 100) % 100;
2620    if (!numberRootThreads) {
2621      if (numberThreads_ < 2)
2622        numberRootThreads = numberModels;
2623      else
2624        numberRootThreads = CoinMin(numberThreads_, numberModels);
2625    }
2626#endif
2627    int otherOptions = (multipleRootTries_ / 10000) % 100;
2628    rootModels = new CbcModel *[numberModels];
2629    unsigned int newSeed = randomSeed_;
2630    if (newSeed == 0) {
2631      double time = fabs(CoinGetTimeOfDay());
2632      while (time >= COIN_INT_MAX)
2633        time *= 0.5;
2634      newSeed = static_cast< unsigned int >(time);
2635    } else if (newSeed < 0) {
2636      newSeed = 123456789;
2637#ifdef COIN_HAS_CLP
2638      OsiClpSolverInterface *clpSolver
2639        = dynamic_cast< OsiClpSolverInterface * >(solver_);
2640      if (clpSolver) {
2641        newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2642      }
2643#endif
2644    }
2645    CoinWarmStartBasis *basis = dynamic_cast< CoinWarmStartBasis * >(solver_->getEmptyWarmStart());
2646    for (int i = 0; i < numberModels; i++) {
2647      rootModels[i] = new CbcModel(*this);
2648      rootModels[i]->setNumberThreads(0);
2649      rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2650      rootModels[i]->setRandomSeed(newSeed + 10000000 * i);
2651      rootModels[i]->randomNumberGenerator()->setSeed(newSeed + 50000000 * i);
2652      rootModels[i]->setMultipleRootTries(0);
2653#ifdef COIN_HAS_NTY
2654      rootModels[i]->zapSymmetry();
2655      rootModels[i]->moreSpecialOptions2_ &= ~(128 | 256); // off nauty
2656#endif
2657      // use seed
2658      rootModels[i]->setSpecialOptions(specialOptions_ | (4194304 | 8388608));
2659      rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & (~(134217728 | 4194304)));
2660      rootModels[i]->setMoreSpecialOptions2(moreSpecialOptions2_ & (~(128 | 256)));
2661      rootModels[i]->solver_->setWarmStart(basis);
2662#ifdef COIN_HAS_CLP
2663      OsiClpSolverInterface *clpSolver
2664        = dynamic_cast< OsiClpSolverInterface * >(rootModels[i]->solver_);
2665#define NEW_RANDOM_BASIS
2666#ifdef NEW_RANDOM_BASIS
2667      if (i == 0)
2668        continue;
2669#endif
2670      if (clpSolver) {
2671        ClpSimplex *simplex = clpSolver->getModelPtr();
2672        if (defaultHandler_)
2673          simplex->setDefaultMessageHandler();
2674        simplex->setRandomSeed(newSeed + 20000000 * i);
2675        simplex->allSlackBasis();
2676        int logLevel = simplex->logLevel();
2677        if (logLevel == 1)
2678          simplex->setLogLevel(0);
2679        if (i != 0) {
2680#ifdef NEW_RANDOM_BASIS
2681          int numberRows = simplex->numberRows();
2682          int throwOut = 20; //2+numberRows/100;
2683          for (int iThrow = 0; iThrow < throwOut; iThrow++) {
2684            double random = simplex->randomNumberGenerator()->randomDouble();
2685            int iStart = static_cast< int >(random * numberRows);
2686            for (int j = iStart; j < numberRows; j++) {
2687              if (simplex->getRowStatus(j) != ClpSimplex::basic) {
2688                simplex->setRowStatus(j, ClpSimplex::basic);
2689                break;
2690              }
2691            }
2692          }
2693          clpSolver->setWarmStart(NULL);
2694#else
2695            double random = simplex->randomNumberGenerator()->randomDouble();
2696            int bias = static_cast< int >(random * (numberIterations / 4));
2697            simplex->setMaximumIterations(numberIterations / 2 + bias);
2698            simplex->primal();
2699            simplex->setMaximumIterations(COIN_INT_MAX);
2700            simplex->dual();
2701#endif
2702        } else {
2703#ifndef NEW_RANDOM_BASIS
2704          simplex->primal();
2705#endif
2706#endif
2707        }
2708#ifdef NEW_RANDOM_BASIS
2709        simplex->setLogLevel(logLevel);
2710        clpSolver->setWarmStart(NULL);
2711#endif
2712      }
2713      for (int j = 0; j < numberHeuristics_; j++)
2714        rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed() + 100000000 * i);
2715      for (int j = 0; j < numberCutGenerators_; j++)
2716        rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
2717    }
2718    delete basis;
2719#ifdef CBC_THREAD
2720    if (numberRootThreads == 1) {
2721#endif
2722      for (int iModel = 0; iModel < numberModels; iModel++) {
2723        doRootCbcThread(rootModels[iModel]);
2724        // see if solved at root node
2725        if (rootModels[iModel]->getMaximumNodes()) {
2726          feasible = false;
2727          break;
2728        }
2729      }
2730#ifdef CBC_THREAD
2731    } else {
2732      Coin_pthread_t *threadId = new Coin_pthread_t[numberRootThreads];
2733      for (int kModel = 0; kModel < numberModels; kModel += numberRootThreads) {
2734        bool finished = false;
2735        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
2736          pthread_create(&(threadId[iModel - kModel].thr), NULL,
2737            doRootCbcThread,
2738            rootModels[iModel]);
2739        }
2740        // wait
2741        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
2742          pthread_join(threadId[iModel - kModel].thr, NULL);
2743        }
2744        // see if solved at root node
2745        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
2746          if (rootModels[iModel]->getMaximumNodes())
2747            finished = true;
2748        }
2749        if (finished) {
2750          feasible = false;
2751          break;
2752        }
2753      }
2754      delete[] threadId;
2755    }
2756#endif
2757    // sort solutions
2758    int *which = new int[numberModels];
2759    double *value = new double[numberModels];
2760    int numberSolutions = 0;
2761    for (int iModel = 0; iModel < numberModels; iModel++) {
2762      if (rootModels[iModel]->bestSolution()) {
2763        which[numberSolutions] = iModel;
2764        value[numberSolutions++] = -rootModels[iModel]->getMinimizationObjValue();
2765      }
2766    }
2767    char general[100];
2768    rootTimeCpu = CoinCpuTime() - rootTimeCpu;
2769    if (numberRootThreads == 1)
2770      sprintf(general, "Multiple root solvers took a total of %.2f seconds\n",
2771        rootTimeCpu);
2772    else
2773      sprintf(general, "Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2774        rootTimeCpu, CoinGetTimeOfDay() - startTimeRoot);
2775    messageHandler()->message(CBC_GENERAL,
2776      messages())
2777      << general << CoinMessageEol;
2778    CoinSort_2(value, value + numberSolutions, which);
2779    // to get name
2780    CbcHeuristicRINS dummyHeuristic;
2781    dummyHeuristic.setHeuristicName("Multiple root solvers");
2782    lastHeuristic_ = &dummyHeuristic;
2783    for (int i = 0; i < numberSolutions; i++) {
2784      double objValue = -value[i];
2785      if (objValue < getCutoff()) {
2786        int iModel = which[i];
2787        setBestSolution(CBC_ROUNDING, objValue,
2788          rootModels[iModel]->bestSolution());
2789      }
2790    }
2791    lastHeuristic_ = NULL;
2792    delete[] which;
2793    delete[] value;
2794  }
2795  // Do heuristics
2796  if (numberObjects_ && !rootModels)
2797    doHeuristicsAtRoot();
2798  if (solverCharacteristics_->solutionAddsCuts()) {
2799    // With some heuristics solver needs a resolve here
2800    solver_->resolve();
2801    if (!isProvenOptimal()) {
2802      solver_->initialSolve();
2803    }
2804  }
2805  /*
2806      Grepping through the code, it would appear that this is a command line
2807      debugging hook.  There's no obvious place in the code where this is set to
2808      a negative value.
2809
2810      User hook, says John.
2811    */
2812  if (intParam_[CbcMaxNumNode] < 0
2813    || numberSolutions_ >= getMaximumSolutions())
2814    eventHappened_ = true; // stop as fast as possible
2815  stoppedOnGap_ = false;
2816  // See if can stop on gap
2817  bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2818  if (canStopOnGap()) {
2819    if (bestPossibleObjective_ < getCutoff())
2820      stoppedOnGap_ = true;
2821    feasible = false;
2822    //eventHappened_=true; // stop as fast as possible
2823  }
2824  /*
2825      Set up for statistics collection, if requested. Standard values are
2826      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2827      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2828      command line debugging hook.
2829    */
2830  statistics_ = NULL;
2831  // Do on switch
2832  if (doStatistics > 0 && doStatistics <= 100) {
2833    maximumStatistics_ = 10000;
2834    statistics_ = new CbcStatistics *[maximumStatistics_];
2835    memset(statistics_, 0, maximumStatistics_ * sizeof(CbcStatistics *));
2836  }
2837  // See if we can add integers
2838  if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_ & 65536) != 0 && !parentModel_ && false) {
2839    int numberIntegers1 = 0;
2840    int numberColumns = solver_->getNumCols();
2841    for (int i = 0; i < numberColumns; i++) {
2842      if (solver_->isInteger(i))
2843        numberIntegers1++;
2844    }
2845    AddIntegers();
2846    // make sure in sync
2847    int numberIntegers2 = 0;
2848    for (int i = 0; i < numberColumns; i++) {
2849      if (solver_->isInteger(i))
2850        numberIntegers2++;
2851    }
2852    if (numberIntegers1 < numberIntegers2) {
2853      findIntegers(true, 2);
2854      convertToDynamic();
2855    }
2856  }
2857
2858  /*
2859      Do an initial round of cut generation for the root node. Depending on the
2860      type of underlying solver, we may want to do this even if the initial query
2861      to the objects indicates they're satisfied.
2862
2863      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2864      loop (including reduced cost fixing) until no cuts are generated, the
2865      change in objective falls off,  or the limit on the number of rounds of cut
2866      generation is exceeded.
2867
2868      At the end of all this, any cuts will be recorded in cuts and also
2869      installed in the solver's constraint system. We'll have reoptimised, and
2870      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2871      adjusted accordingly).
2872    */
2873  int iObject;
2874  int numberUnsatisfied = 0;
2875  delete[] currentSolution_;
2876  currentSolution_ = new double[numberColumns];
2877  testSolution_ = currentSolution_;
2878  memcpy(currentSolution_, solver_->getColSolution(),
2879    numberColumns * sizeof(double));
2880  // point to useful information
2881  OsiBranchingInformation usefulInfo = usefulInformation();
2882
2883  for (iObject = 0; iObject < numberObjects_; iObject++) {
2884    double infeasibility = object_[iObject]->checkInfeasibility(&usefulInfo);
2885    if (infeasibility)
2886      numberUnsatisfied++;
2887  }
2888  // replace solverType
2889  double *tightBounds = NULL;
2890  if (solverCharacteristics_->tryCuts()) {
2891
2892    if (numberUnsatisfied) {
2893      // User event
2894      if (!eventHappened_ && feasible) {
2895        if (rootModels) {
2896          // for fixings
2897          int numberColumns = solver_->getNumCols();
2898          tightBounds = new double[2 * numberColumns];
2899          {
2900            const double *lower = solver_->getColLower();
2901            const double *upper = solver_->getColUpper();
2902            for (int i = 0; i < numberColumns; i++) {
2903              tightBounds[2 * i + 0] = lower[i];
2904              tightBounds[2 * i + 1] = upper[i];
2905            }
2906          }
2907          int numberModels = multipleRootTries_ % 100;
2908          const OsiSolverInterface **solvers = new const OsiSolverInterface *[numberModels];
2909          int numberRows = continuousSolver_->getNumRows();
2910          int maxCuts = 0;
2911          for (int i = 0; i < numberModels; i++) {
2912            solvers[i] = rootModels[i]->solver();
2913            const double *lower = solvers[i]->getColLower();
2914            const double *upper = solvers[i]->getColUpper();
2915            for (int j = 0; j < numberColumns; j++) {
2916              tightBounds[2 * j + 0] = CoinMax(lower[j], tightBounds[2 * j + 0]);
2917              tightBounds[2 * j + 1] = CoinMin(upper[j], tightBounds[2 * j + 1]);
2918            }
2919            int numberRows2 = solvers[i]->getNumRows();
2920            assert(numberRows2 >= numberRows);
2921            maxCuts += numberRows2 - numberRows;
2922            // accumulate statistics
2923            for (int j = 0; j < numberCutGenerators_; j++) {
2924              generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2925            }
2926          }
2927          for (int j = 0; j < numberCutGenerators_; j++) {
2928            generator_[j]->scaleBackStatistics(numberModels);
2929          }
2930          //CbcRowCuts rowCut(maxCuts);
2931          const OsiRowCutDebugger *debugger = NULL;
2932          if ((specialOptions_ & 1) != 0)
2933            debugger = solver_->getRowCutDebugger();
2934          for (int iModel = 0; iModel < numberModels; iModel++) {
2935            int numberRows2 = solvers[iModel]->getNumRows();
2936            const CoinPackedMatrix *rowCopy = solvers[iModel]->getMatrixByRow();
2937            const int *rowLength = rowCopy->getVectorLengths();
2938            const double *elements = rowCopy->getElements();
2939            const int *column = rowCopy->getIndices();
2940            const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
2941            const double *rowLower = solvers[iModel]->getRowLower();
2942            const double *rowUpper = solvers[iModel]->getRowUpper();
2943            for (int iRow = numberRows; iRow < numberRows2; iRow++) {
2944              OsiRowCut rc;
2945              rc.setLb(rowLower[iRow]);
2946              rc.setUb(rowUpper[iRow]);
2947              CoinBigIndex start = rowStart[iRow];
2948              rc.setRow(rowLength[iRow], column + start, elements + start, false);
2949              if (debugger)
2950                CoinAssert(!debugger->invalidCut(rc));
2951              globalCuts_.addCutIfNotDuplicate(rc);
2952            }
2953            //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
2954            //numberCuts += cutsAdded;
2955            //printf("Model %d gave %d cuts (out of %d possible)\n",
2956            //     iModel,cutsAdded,numberRows2-numberRows);
2957          }
2958          // normally replace global cuts
2959          //if (!globalCuts_.())
2960          //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2961          //rowCut.addCuts(globalCuts_);
2962          int nTightened = 0;
2963          assert(feasible);
2964          {
2965            double tolerance = 1.0e-5;
2966            const double *lower = solver_->getColLower();
2967            const double *upper = solver_->getColUpper();
2968            for (int i = 0; i < numberColumns; i++) {
2969              if (tightBounds[2 * i + 0] > tightBounds[2 * i + 1] + 1.0e-9) {
2970                feasible = false;
2971                char general[200];
2972                sprintf(general, "Solvers give infeasible bounds on %d %g,%g was %g,%g - search finished\n",
2973                  i, tightBounds[2 * i + 0], tightBounds[2 * i + 1], lower[i], upper[i]);
2974                messageHandler()->message(CBC_GENERAL, messages())
2975                  << general << CoinMessageEol;
2976                break;
2977              }
2978              double oldLower = lower[i];
2979              double oldUpper = upper[i];
2980              if (tightBounds[2 * i + 0] > oldLower + tolerance) {
2981                nTightened++;
2982                solver_->setColLower(i, tightBounds[2 * i + 0]);
2983              }
2984              if (tightBounds[2 * i + 1] < oldUpper - tolerance) {
2985                nTightened++;
2986                solver_->setColUpper(i, tightBounds[2 * i + 1]);
2987              }
2988            }
2989          }
2990          delete[] tightBounds;
2991          tightBounds = NULL;
2992          char printBuffer[200];
2993          sprintf(printBuffer, "%d solvers added %d different cuts out of pool of %d",
2994            numberModels, globalCuts_.sizeRowCuts(), maxCuts);
2995          messageHandler()->message(CBC_GENERAL, messages())
2996            << printBuffer << CoinMessageEol;
2997          if (nTightened) {
2998            sprintf(printBuffer, "%d bounds were tightened",
2999              nTightened);
3000            messageHandler()->message(CBC_GENERAL, messages())
3001              << printBuffer << CoinMessageEol;
3002          }
3003          delete[] solvers;
3004        }
3005        if (!parentModel_ && (moreSpecialOptions_ & 67108864) != 0) {
3006          // load cuts from file
3007          FILE *fp = fopen("global.cuts", "rb");
3008          if (fp) {
3009            size_t nRead;
3010            int numberColumns = solver_->getNumCols();
3011            int numCols;
3012            nRead = fread(&numCols, sizeof(int), 1, fp);
3013            if (nRead != 1)
3014              throw("Error in fread");
3015            if (numberColumns != numCols) {
3016              printf("Mismatch on columns %d %d\n", numberColumns, numCols);
3017              fclose(fp);
3018            } else {
3019              // If rootModel just do some
3020              double threshold = -1.0;
3021              if (!multipleRootTries_)
3022                threshold = 0.5;
3023              int initialCuts = 0;
3024              int initialGlobal = globalCuts_.sizeRowCuts();
3025              double *elements = new double[numberColumns + 2];
3026              int *indices = new int[numberColumns];
3027              int numberEntries = 1;
3028              while (numberEntries > 0) {
3029                nRead = fread(&numberEntries, sizeof(int), 1, fp);
3030                if (nRead != 1)
3031                  throw("Error in fread");
3032                double randomNumber = randomNumberGenerator_.randomDouble();
3033                if (numberEntries > 0) {
3034                  initialCuts++;
3035                  nRead = fread(elements, sizeof(double), numberEntries + 2, fp);
3036                  if (nRead != static_cast< size_t >(numberEntries + 2))
3037                    throw("Error in fread");
3038                  nRead = fread(indices, sizeof(int), numberEntries, fp);
3039                  if (nRead != static_cast< size_t >(numberEntries))
3040                    throw("Error in fread");
3041                  if (randomNumber > threshold) {
3042                    OsiRowCut rc;
3043                    rc.setLb(elements[numberEntries]);
3044                    rc.setUb(elements[numberEntries + 1]);
3045                    rc.setRow(numberEntries, indices, elements,
3046                      false);
3047                    rc.setGloballyValidAsInteger(2);
3048                    globalCuts_.addCutIfNotDuplicate(rc);
3049                  }
3050                }
3051              }
3052              fclose(fp);
3053              // fixes
3054              int nTightened = 0;
3055              fp = fopen("global.fix", "rb");
3056              if (fp) {
3057                nRead = fread(indices, sizeof(int), 2, fp);
3058                if (nRead != 2)
3059                  throw("Error in fread");
3060                if (numberColumns != indices[0]) {
3061                  printf("Mismatch on columns %d %d\n", numberColumns,
3062                    indices[0]);
3063                } else {
3064                  indices[0] = 1;
3065                  while (indices[0] >= 0) {
3066                    nRead = fread(indices, sizeof(int), 2, fp);
3067                    if (nRead != 2)
3068                      throw("Error in fread");
3069                    int iColumn = indices[0];
3070                    if (iColumn >= 0) {
3071                      nTightened++;
3072                      nRead = fread(elements, sizeof(double), 4, fp);
3073                      if (nRead != 4)
3074                        throw("Error in fread");
3075                      solver_->setColLower(iColumn, elements[0]);
3076                      solver_->setColUpper(iColumn, elements[1]);
3077                    }
3078                  }
3079                }
3080              }
3081              if (fp)
3082                fclose(fp);
3083              char printBuffer[200];
3084              sprintf(printBuffer, "%d cuts read in of which %d were unique, %d bounds tightened",
3085                initialCuts,
3086                globalCuts_.sizeRowCuts() - initialGlobal, nTightened);
3087              messageHandler()->message(CBC_GENERAL, messages())
3088                << printBuffer << CoinMessageEol;
3089              delete[] elements;
3090              delete[] indices;
3091            }
3092          }
3093        }
3094        if (feasible)
3095          feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3096            NULL);
3097        if (multipleRootTries_ && (moreSpecialOptions_ & 134217728) != 0) {
3098          FILE *fp = NULL;
3099          size_t nRead;
3100          int numberColumns = solver_->getNumCols();
3101          int initialCuts = 0;
3102          if ((moreSpecialOptions_ & 134217728) != 0) {
3103            // append so go down to end
3104            fp = fopen("global.cuts", "r+b");
3105            if (fp) {
3106              int numCols;
3107              nRead = fread(&numCols, sizeof(int), 1, fp);
3108              if (nRead != 1)
3109                throw("Error in fread");
3110              if (numberColumns != numCols) {
3111                printf("Mismatch on columns %d %d\n", numberColumns, numCols);
3112                fclose(fp);
3113                fp = NULL;
3114              }
3115            }
3116          }
3117          double *elements = new double[numberColumns + 2];
3118          int *indices = new int[numberColumns];
3119          if (fp) {
3120            int numberEntries = 1;
3121            while (numberEntries > 0) {
3122              fpos_t position;
3123              fgetpos(fp, &position);
3124              nRead = fread(&numberEntries, sizeof(int), 1, fp);
3125              if (nRead != 1)
3126                throw("Error in fread");
3127              if (numberEntries > 0) {
3128                initialCuts++;
3129                nRead = fread(elements, sizeof(double), numberEntries + 2, fp);
3130                if (nRead != static_cast< size_t >(numberEntries + 2))
3131                  throw("Error in fread");
3132                nRead = fread(indices, sizeof(int), numberEntries, fp);
3133                if (nRead != static_cast< size_t >(numberEntries))
3134                  throw("Error in fread");
3135              } else {
3136                // end
3137                fsetpos(fp, &position);
3138              }
3139            }
3140          } else {
3141            fp = fopen("global.cuts", "wb");
3142            size_t nWrite;
3143            nWrite = fwrite(&numberColumns, sizeof(int), 1, fp);
3144            if (nWrite != 1)
3145              throw("Error in fwrite");
3146          }
3147          size_t nWrite;
3148          // now append binding cuts
3149          int numberC = continuousSolver_->getNumRows();
3150          int numberRows = solver_->getNumRows();
3151          printf("Saving %d cuts (up from %d)\n",
3152            initialCuts + numberRows - numberC, initialCuts);
3153          const double *rowLower = solver_->getRowLower();
3154          const double *rowUpper = solver_->getRowUpper();
3155          // Row copy
3156          CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3157          const double *elementByRow = matrixByRow.getElements();
3158          const int *column = matrixByRow.getIndices();
3159          const CoinBigIndex *rowStart = matrixByRow.getVectorStarts();
3160          const int *rowLength = matrixByRow.getVectorLengths();
3161          for (int iRow = numberC; iRow < numberRows; iRow++) {
3162            int n = rowLength[iRow];
3163            assert(n);
3164            CoinBigIndex start = rowStart[iRow];
3165            memcpy(elements, elementByRow + start, n * sizeof(double));
3166            memcpy(indices, column + start, n * sizeof(int));
3167            elements[n] = rowLower[iRow];
3168            elements[n + 1] = rowUpper[iRow];
3169            nWrite = fwrite(&n, sizeof(int), 1, fp);
3170            if (nWrite != 1)
3171              throw("Error in fwrite");
3172            nWrite = fwrite(elements, sizeof(double), n + 2, fp);
3173            if (nWrite != static_cast< size_t >(n + 2))
3174              throw("Error in fwrite");
3175            nWrite = fwrite(indices, sizeof(int), n, fp);
3176            if (nWrite != static_cast< size_t >(n))
3177              throw("Error in fwrite");
3178          }
3179          // eof marker
3180          int eofMarker = -1;
3181          nWrite = fwrite(&eofMarker, sizeof(int), 1, fp);
3182          if (nWrite != 1)
3183            throw("Error in fwrite");
3184          fclose(fp);
3185          // do tighter bounds (? later extra to original columns)
3186          int nTightened = 0;
3187          const double *lower = solver_->getColLower();
3188          const double *upper = solver_->getColUpper();
3189          const double *originalLower = continuousSolver_->getColLower();
3190          const double *originalUpper = continuousSolver_->getColUpper();
3191          double tolerance = 1.0e-5;
3192          for (int i = 0; i < numberColumns; i++) {
3193            if (lower[i] > originalLower[i] + tolerance) {
3194              nTightened++;
3195            }
3196            if (upper[i] < originalUpper[i] - tolerance) {
3197              nTightened++;
3198            }
3199          }
3200          if (nTightened) {
3201            fp = fopen("global.fix", "wb");
3202            size_t nWrite;
3203            indices[0] = numberColumns;
3204            if (originalColumns_)
3205              indices[1] = COIN_INT_MAX;
3206            else
3207              indices[1] = -1;
3208            nWrite = fwrite(indices, sizeof(int), 2, fp);
3209            if (nWrite != 2)
3210              throw("Error in fwrite");
3211            for (int i = 0; i < numberColumns; i++) {
3212              int nTightened = 0;
3213              if (lower[i] > originalLower[i] + tolerance) {
3214                nTightened++;
3215              }
3216              if (upper[i] < originalUpper[i] - tolerance) {
3217                nTightened++;
3218              }
3219              if (nTightened) {
3220                indices[0] = i;
3221                if (originalColumns_)
3222                  indices[1] = originalColumns_[i];
3223                elements[0] = lower[i];
3224                elements[1] = upper[i];
3225                elements[2] = originalLower[i];
3226                elements[3] = originalUpper[i];
3227                nWrite = fwrite(indices, sizeof(int), 2, fp);
3228                if (nWrite != 2)
3229                  throw("Error in fwrite");
3230                nWrite = fwrite(elements, sizeof(double), 4, fp);
3231                if (nWrite != 4)
3232                  throw("Error in fwrite");
3233              }
3234            }
3235            // eof marker
3236            indices[0] = -1;
3237            nWrite = fwrite(indices, sizeof(int), 2, fp);
3238            if (nWrite != 2)
3239              throw("Error in fwrite");
3240            fclose(fp);
3241          }
3242          delete[] elements;
3243          delete[] indices;
3244        }
3245        if ((specialOptions_ & 524288) != 0 && !parentModel_
3246          && storedRowCuts_) {
3247          if (feasible) {
3248            /* pick up stuff and try again
3249                        add cuts, maybe keep around
3250                        do best solution and if so new heuristics
3251                        obviously tighten bounds
3252                        */
3253            // A and B probably done on entry
3254            // A) tight bounds on integer variables
3255            const double *lower = solver_->getColLower();
3256            const double *upper = solver_->getColUpper();
3257            const double *tightLower = storedRowCuts_->tightLower();
3258            const double *tightUpper = storedRowCuts_->tightUpper();
3259            int nTightened = 0;
3260            for (int i = 0; i < numberIntegers_; i++) {
3261              int iColumn = integerVariable_[i];
3262              if (tightLower[iColumn] > lower[iColumn]) {
3263                nTightened++;
3264                solver_->setColLower(iColumn, tightLower[iColumn]);
3265              }
3266              if (tightUpper[iColumn] < upper[iColumn]) {
3267                nTightened++;
3268                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3269              }
3270            }
3271            if (nTightened)
3272              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3273            if (storedRowCuts_->bestObjective() < bestObjective_) {
3274              // B) best solution
3275              double objValue = storedRowCuts_->bestObjective();
3276              setBestSolution(CBC_SOLUTION, objValue,
3277                storedRowCuts_->bestSolution());
3278              // Do heuristics
3279              // Allow RINS
3280              for (int i = 0; i < numberHeuristics_; i++) {
3281                CbcHeuristicRINS *rins
3282                  = dynamic_cast< CbcHeuristicRINS * >(heuristic_[i]);
3283                if (rins) {
3284                  rins->setLastNode(-100);
3285                }
3286              }
3287              doHeuristicsAtRoot();
3288            }
3289#ifdef JJF_ZERO
3290            int nCuts = storedRowCuts_->sizeRowCuts();
3291            // add to global list
3292            for (int i = 0; i < nCuts; i++) {
3293              OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3294              newCut.setGloballyValidAsInteger(2);
3295              newCut.mutableRow().setTestForDuplicateIndex(false);
3296              globalCuts_.insert(newCut);
3297            }
3298#else
3299          addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3300            true, false, false, -200);
3301#endif
3302            // Set cuts as active
3303            delete[] addedCuts_;
3304            maximumNumberCuts_ = cuts.sizeRowCuts();
3305            if (maximumNumberCuts_) {
3306              addedCuts_ = new CbcCountRowCut *[maximumNumberCuts_];
3307            } else {
3308              addedCuts_ = NULL;
3309            }
3310            for (int i = 0; i < maximumNumberCuts_; i++)
3311              addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3312                NULL, -1, -1, 2);
3313            COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3314            cuts = OsiCuts();
3315            currentNumberCuts_ = maximumNumberCuts_;
3316            feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3317              NULL);
3318            for (int i = 0; i < maximumNumberCuts_; i++)
3319              delete addedCuts_[i];
3320          }
3321          delete storedRowCuts_;
3322          storedRowCuts_ = NULL;
3323        }
3324      } else {
3325        feasible = false;
3326      }
3327    } else if (solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3328      // may generate cuts and turn the solution
3329      //to an infeasible one
3330      feasible = solveWithCuts(cuts, 2,
3331        NULL);
3332    }
3333  }
3334  if (rootModels) {
3335    int numberModels = multipleRootTries_ % 100;
3336    for (int i = 0; i < numberModels; i++)
3337      delete rootModels[i];
3338    delete[] rootModels;
3339  }
3340  // check extra info on feasibility
3341  if (!solverCharacteristics_->mipFeasible())
3342    feasible = false;
3343  // If max nodes==0 - don't do strong branching
3344  if (!getMaximumNodes()) {
3345    if (feasible)
3346      feasible = false;
3347    else
3348      setMaximumNodes(1); //allow to stop on success
3349  }
3350  topOfTree_ = NULL;
3351#ifdef CLP_RESOLVE
3352  if ((moreSpecialOptions_ & 2097152) != 0 && !parentModel_ && feasible) {
3353    OsiClpSolverInterface *clpSolver
3354      = dynamic_cast< OsiClpSolverInterface * >(solver_);
3355    if (clpSolver)
3356      resolveClp(clpSolver, 0);
3357  }
3358#endif
3359  // make cut generators less aggressive
3360  for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3361    CglCutGenerator *generator = generator_[iCutGenerator]->generator();
3362    generator->setAggressiveness(generator->getAggressiveness() - 100);
3363  }
3364  currentNumberCuts_ = numberNewCuts_;
3365  if (solverCharacteristics_->solutionAddsCuts()) {
3366    // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3367    solver_->resolve();
3368    if (!isProvenOptimal()) {
3369      solver_->initialSolve();
3370    }
3371  }
3372  // See if can stop on gap
3373  bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3374  if (canStopOnGap()) {
3375    if (bestPossibleObjective_ < getCutoff())
3376      stoppedOnGap_ = true;
3377    feasible = false;
3378  }
3379  // User event
3380  if (eventHappened_)
3381    feasible = false;
3382#if defined(COIN_HAS_CLP) && defined(COIN_HAS_CPX)
3383  /*
3384      This is the notion of using Cbc stuff to get going, then calling cplex to
3385      finish off.
3386    */
3387  if (feasible && (specialOptions_ & 16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3388    // Use Cplex to do search!
3389    double time1 = CoinCpuTime();
3390    OsiClpSolverInterface *clpSolver
3391      = dynamic_cast< OsiClpSolverInterface * >(solver_);
3392    OsiCpxSolverInterface cpxSolver;
3393    double direction = clpSolver->getObjSense();
3394    cpxSolver.setObjSense(direction);
3395    // load up cplex
3396    const CoinPackedMatrix *matrix = continuousSolver_->getMatrixByCol();
3397    const double *rowLower = continuousSolver_->getRowLower();
3398    const double *rowUpper = continuousSolver_->getRowUpper();
3399    const double *columnLower = continuousSolver_->getColLower();
3400    const double *columnUpper = continuousSolver_->getColUpper();
3401    const double *objective = continuousSolver_->getObjCoefficients();
3402    cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3403      objective, rowLower, rowUpper);
3404    double *setSol = new double[numberIntegers_];
3405    int *setVar = new int[numberIntegers_];
3406    // cplex doesn't know about objective offset
3407    double offset = clpSolver->getModelPtr()->objectiveOffset();
3408    for (int i = 0; i < numberIntegers_; i++) {
3409      int iColumn = integerVariable_[i];
3410      cpxSolver.setInteger(iColumn);
3411      if (bestSolution_) {
3412        setSol[i] = bestSolution_[iColumn];
3413        setVar[i] = iColumn;
3414      }
3415    }
3416    CPXENVptr env = cpxSolver.getEnvironmentPtr();
3417    CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3418    cpxSolver.switchToMIP();
3419    if (bestSolution_) {
3420#if 0
3421            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3422#else
3423        int zero = 0;
3424        CPXaddmipstarts(env, lpPtr, 1, numberIntegers_, &zero, setVar, setSol, NULL, NULL);
3425#endif
3426    }
3427    if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3428      // add cuts
3429      const CoinPackedMatrix *matrix = clpSolver->getMatrixByRow();
3430      const double *rhs = clpSolver->getRightHandSide();
3431      const char *rowSense = clpSolver->getRowSense();
3432      const double *elementByRow = matrix->getElements();
3433      const int *column = matrix->getIndices();
3434      const CoinBigIndex *rowStart = matrix->getVectorStarts();
3435      const int *rowLength = matrix->getVectorLengths();
3436      int nStart = continuousSolver_->getNumRows();
3437      int nRows = clpSolver->getNumRows();
3438      int size = rowStart[nRows - 1] + rowLength[nRows - 1] - rowStart[nStart];
3439      int nAdd = 0;
3440      double *rmatval = new double[size];
3441      int *rmatind = new int[size];
3442      int *rmatbeg = new int[nRows - nStart + 1];
3443      size = 0;
3444      rmatbeg[0] = 0;
3445      for (int i = nStart; i < nRows; i++) {
3446        for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3447          rmatind[size] = column[k];
3448          rmatval[size++] = elementByRow[k];
3449        }
3450        nAdd++;
3451        rmatbeg[nAdd] = size;
3452      }
3453      CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3454        rhs, rowSense, rmatbeg,
3455        rmatind, rmatval, NULL);
3456      CPXsetintparam(env, CPX_PARAM_REDUCE,
3457        // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3458        CPX_PREREDUCE_PRIMALONLY);
3459    }
3460    if (getCutoff() < 1.0e50) {
3461      double useCutoff = getCutoff() + offset;
3462      if (bestObjective_ < 1.0e50)
3463        useCutoff = bestObjective_ + offset + 1.0e-7;
3464      cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff * direction);
3465      if (direction > 0.0)
3466        CPXsetdblparam(env, CPX_PARAM_CUTUP, useCutoff); // min
3467      else
3468        CPXsetdblparam(env, CPX_PARAM_CUTLO, useCutoff); // max
3469    }
3470    CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3471    delete[] setSol;
3472    delete[] setVar;
3473    char printBuffer[200];
3474    if (offset) {
3475      sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3476        -offset);
3477      messageHandler()->message(CBC_GENERAL, messages())
3478        << printBuffer << CoinMessageEol;
3479      cpxSolver.setDblParam(OsiObjOffset, offset);
3480    }
3481    cpxSolver.branchAndBound();
3482    double timeTaken = CoinCpuTime() - time1;
3483    sprintf(printBuffer, "Cplex took %g seconds",
3484      timeTaken);
3485    messageHandler()->message(CBC_GENERAL, messages())
3486      << printBuffer << CoinMessageEol;
3487    numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3488    numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3489    double value = cpxSolver.getObjValue() * direction;
3490    if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3491      feasible = true;
3492      clpSolver->setWarmStart(NULL);
3493      // try and do solution
3494      double *newSolution = CoinCopyOfArray(cpxSolver.getColSolution(),
3495        getNumCols());
3496      setBestSolution(CBC_STRONGSOL, value, newSolution);
3497      delete[] newSolution;
3498    }
3499    feasible = false;
3500  }
3501#endif
3502  if (!parentModel_ && (moreSpecialOptions_ & 268435456) != 0) {
3503    // try idiotic idea
3504    CbcObject *obj = new CbcIdiotBranch(this);
3505    obj->setPriority(1); // temp
3506    addObjects(1, &obj);
3507    delete obj;
3508  }
3509
3510  /*
3511      A hook to use clp to quickly explore some part of the tree.
3512    */
3513  if (fastNodeDepth_ == 1000 && /*!parentModel_*/ (specialOptions_ & 2048) == 0) {
3514    fastNodeDepth_ = -1;
3515    CbcObject *obj = new CbcFollowOn(this);
3516    obj->setPriority(1);
3517    addObjects(1, &obj);
3518    delete obj;
3519  }
3520  int saveNumberSolves = numberSolves_;
3521  int saveNumberIterations = numberIterations_;
3522  if ((fastNodeDepth_ >= 0 || (moreSpecialOptions_ & 33554432) != 0)
3523    && /*!parentModel_*/ (specialOptions_ & 2048) == 0) {
3524    // add in a general depth object doClp
3525    int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3526    if ((moreSpecialOptions_ & 33554432) != 0)
3527      type = 12;
3528    else
3529      fastNodeDepth_ += 1000000; // mark as done
3530    CbcObject *obj = new CbcGeneralDepth(this, type);
3531    addObjects(1, &obj);
3532    delete obj;
3533    // fake number of objects
3534    numberObjects_--;
3535    if (parallelMode() < -1) {
3536      // But make sure position is correct
3537      OsiObject *obj2 = object_[numberObjects_];
3538      obj = dynamic_cast< CbcObject * >(obj2);
3539      assert(obj);
3540      obj->setPosition(numberObjects_);
3541    }
3542  }
3543#ifdef COIN_HAS_CLP
3544#ifdef NO_CRUNCH
3545  if (true) {
3546    OsiClpSolverInterface *clpSolver
3547      = dynamic_cast< OsiClpSolverInterface * >(solver_);
3548    if (clpSolver && !parentModel_) {
3549      ClpSimplex *clpSimplex = clpSolver->getModelPtr();
3550      clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
3551      //clpSimplex->startPermanentArrays();
3552      clpSimplex->setPersistenceFlag(2);
3553    }
3554  }
3555#endif
3556#endif
3557  // Save copy of solver
3558  OsiSolverInterface *saveSolver = NULL;
3559  if (!parentModel_ && (specialOptions_ & (512 + 32768)) != 0)
3560    saveSolver = solver_->clone();
3561  double checkCutoffForRestart = 1.0e100;
3562  saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3563  if ((specialOptions_ & 262144) != 0 && !parentModel_) {
3564    // Save stuff and return!
3565    storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3566      solver_->getColLower(),
3567      solver_->getColUpper());
3568    delete[] lowerBefore;
3569    delete[] upperBefore;
3570    delete saveSolver;
3571    if (flipObjective)
3572      flipModel();
3573    return;
3574  }
3575  /*
3576      We've taken the continuous relaxation as far as we can. Time to branch.
3577      The first order of business is to actually create a node. chooseBranch
3578      currently uses strong branching to evaluate branch object candidates,
3579      unless forced back to simple branching. If chooseBranch concludes that a
3580      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3581      == -2) when forced to integer values, it returns here immediately.
3582
3583      Monotone variables trigger a call to resolve(). If the problem remains
3584      feasible, try again to choose a branching variable. At the end of the loop,
3585      resolved == true indicates that some variables were fixed.
3586
3587      Loss of feasibility will result in the deletion of newNode.
3588    */
3589
3590  bool resolved = false;
3591  CbcNode *newNode = NULL;
3592  numberFixedAtRoot_ = 0;
3593  numberFixedNow_ = 0;
3594  if (!parentModel_ && (moreSpecialOptions2_ & 2) != 0) {
3595#ifdef COIN_HAS_CLP
3596    OsiClpSolverInterface *clpSolver
3597      = dynamic_cast< OsiClpSolverInterface * >(solver_);
3598    if (clpSolver) {
3599      if (getCutoff() > 1.0e20) {
3600        printf("Zapping costs\n");
3601        int numberColumns = solver_->getNumCols();
3602        double *zeroCost = new double[numberColumns];
3603        // could make small random
3604        memset(zeroCost, 0, numberColumns * sizeof(double));
3605        solver_->setObjective(zeroCost);
3606        double objValue = solver_->getObjValue();
3607        solver_->setDblParam(OsiObjOffset, -objValue);
3608        clpSolver->getModelPtr()->setObjectiveValue(objValue);
3609        delete[] zeroCost;
3610      } else {
3611        moreSpecialOptions2_ &= ~2;
3612      }
3613    } else {
3614#endif
3615      moreSpecialOptions2_ &= ~2;
3616#ifdef COIN_HAS_CLP
3617    }
3618#endif
3619  }
3620  int numberIterationsAtContinuous = numberIterations_;
3621  //solverCharacteristics_->setSolver(solver_);
3622  if (feasible) {
3623    // mark all cuts as globally valid
3624    int numberCuts = cuts.sizeRowCuts();
3625    resizeWhichGenerator(0, numberCuts);
3626    for (int i = 0; i < numberCuts; i++) {
3627      cuts.rowCutPtr(i)->setGloballyValid();
3628      whichGenerator_[i] = 20000 + (whichGenerator_[i] % 10000);
3629    }
3630#define HOTSTART -1
3631#if HOTSTART < 0
3632    if (bestSolution_ && !parentModel_ && !hotstartSolution_ && (moreSpecialOptions_ & 1024) != 0 && (specialOptions_ & 2048) == 0) {
3633      // Set priorities so only branch on ones we need to
3634      // use djs and see if only few branches needed
3635#ifndef NDEBUG
3636      double integerTolerance = getIntegerTolerance();
3637#endif
3638      bool possible = true;
3639      const double *saveLower = continuousSolver_->getColLower();
3640      const double *saveUpper = continuousSolver_->getColUpper();
3641      for (int i = 0; i < numberObjects_; i++) {
3642        const CbcSimpleInteger *thisOne = dynamic_cast< const CbcSimpleInteger * >(object_[i]);
3643        if (thisOne) {
3644          int iColumn = thisOne->columnNumber();
3645          if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3646            possible = false;
3647            break;
3648          }
3649        } else {
3650          possible = false;
3651          break;
3652        }
3653      }
3654      if (possible) {
3655        OsiSolverInterface *solver = continuousSolver_->clone();
3656        int numberColumns = solver->getNumCols();
3657        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3658          double value = bestSolution_[iColumn];
3659          value = CoinMax(value, saveLower[iColumn]);
3660          value = CoinMin(value, saveUpper[iColumn]);
3661          value = floor(value + 0.5);
3662          if (solver->isInteger(iColumn)) {
3663            solver->setColLower(iColumn, value);
3664            solver->setColUpper(iColumn, value);
3665          }
3666        }
3667        solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3668        // objlim and all slack
3669        double direction = solver->getObjSense();
3670        solver->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
3671        CoinWarmStartBasis *basis = dynamic_cast< CoinWarmStartBasis * >(solver->getEmptyWarmStart());
3672        solver->setWarmStart(basis);
3673        delete basis;
3674        bool changed = true;
3675        hotstartPriorities_ = new int[numberColumns];
3676        for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3677          hotstartPriorities_[iColumn] = 1;
3678        while (changed) {
3679          changed = false;
3680          solver->resolve();
3681          if (!solver->isProvenOptimal()) {
3682            possible = false;
3683            break;
3684          }
3685          const double *dj = solver->getReducedCost();
3686          const double *colLower = solver->getColLower();
3687          const double *colUpper = solver->getColUpper();
3688          const double *solution = solver->getColSolution();
3689          int nAtLbNatural = 0;
3690          int nAtUbNatural = 0;
3691          int nZeroDj = 0;
3692          int nForced = 0;
3693          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3694            double value = solution[iColumn];
3695            value = CoinMax(value, saveLower[iColumn]);
3696            value = CoinMin(value, saveUpper[iColumn]);
3697            if (solver->isInteger(iColumn)) {
3698              assert(fabs(value - solution[iColumn]) <= integerTolerance);
3699              if (hotstartPriorities_[iColumn] == 1) {
3700                if (dj[iColumn] < -1.0e-6) {
3701                  // negative dj
3702                  if (saveUpper[iColumn] == colUpper[iColumn]) {
3703                    nAtUbNatural++;
3704                    hotstartPriorities_[iColumn] = 2;
3705                    solver->setColLower(iColumn, saveLower[iColumn]);
3706                    solver->setColUpper(iColumn, saveUpper[iColumn]);
3707                  } else {
3708                    nForced++;
3709                  }
3710                } else if (dj[iColumn] > 1.0e-6) {
3711                  // positive dj
3712                  if (saveLower[iColumn] == colLower[iColumn]) {
3713                    nAtLbNatural++;
3714                    hotstartPriorities_[iColumn] = 2;
3715                    solver->setColLower(iColumn, saveLower[iColumn]);
3716                    solver->setColUpper(iColumn, saveUpper[iColumn]);
3717                  } else {
3718                    nForced++;
3719                  }
3720                } else {
3721                  // zero dj
3722                  nZeroDj++;
3723                }
3724              }
3725            }
3726          }
3727#if CBC_USEFUL_PRINTING > 1
3728          printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3729            nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3730#endif
3731          if (nAtLbNatural || nAtUbNatural) {
3732            changed = true;
3733          } else {
3734            if (nForced + nZeroDj > 5000 || (nForced + nZeroDj) * 2 > numberIntegers_)
3735              possible = false;
3736          }
3737        }
3738        delete solver;
3739      }
3740      if (possible) {
3741        setHotstartSolution(bestSolution_);
3742        if (!saveCompare) {
3743          // create depth first comparison
3744          saveCompare = nodeCompare_;
3745          // depth first
3746          nodeCompare_ = new CbcCompareDepth();
3747          tree_->setComparison(*nodeCompare_);
3748        }
3749      } else {
3750        delete[] hotstartPriorities_;
3751        hotstartPriorities_ = NULL;
3752      }
3753    }
3754#endif
3755#if HOTSTART > 0
3756    if (hotstartSolution_ && !hotstartPriorities_) {
3757      // Set up hot start
3758      OsiSolverInterface *solver = solver_->clone();
3759      double direction = solver_->getObjSense();
3760      int numberColumns = solver->getNumCols();
3761      double *saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3762      double *saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3763      // move solution
3764      solver->setColSolution(hotstartSolution_);
3765      // point to useful information
3766      const double *saveSolution = testSolution_;
3767      testSolution_ = solver->getColSolution();
3768      OsiBranchingInformation usefulInfo = usefulInformation();
3769      testSolution_ = saveSolution;
3770      /*
3771            Run through the objects and use feasibleRegion() to set variable bounds
3772            so as to fix the variables specified in the objects at their value in this
3773            solution. Since the object list contains (at least) one object for every
3774            integer variable, this has the effect of fixing all integer variables.
3775            */
3776      for (int i = 0; i < numberObjects_; i++)
3777        object_[i]->feasibleRegion(solver, &usefulInfo);
3778      solver->resolve();
3779      assert(solver->isProvenOptimal());
3780      double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0);
3781      const double *dj = solver->getReducedCost();
3782      const double *colLower = solver->getColLower();
3783      const double *colUpper = solver->getColUpper();
3784      const double *solution = solver->getColSolution();
3785      int nAtLbNatural = 0;
3786      int nAtUbNatural = 0;
3787      int nAtLbNaturalZero = 0;
3788      int nAtUbNaturalZero = 0;
3789      int nAtLbFixed = 0;
3790      int nAtUbFixed = 0;
3791      int nAtOther = 0;
3792      int nAtOtherNatural = 0;
3793      int nNotNeeded = 0;
3794      delete[] hotstartSolution_;
3795      hotstartSolution_ = new double[numberColumns];
3796      delete[] hotstartPriorities_;
3797      hotstartPriorities_ = new int[numberColumns];
3798      int *order = (int *)saveUpper;
3799      int nFix = 0;
3800      double bestRatio = COIN_DBL_MAX;
3801      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3802        double value = solution[iColumn];
3803        value = CoinMax(value, saveLower[iColumn]);
3804        value = CoinMin(value, saveUpper[iColumn]);
3805        double sortValue = COIN_DBL_MAX;
3806        if (solver->isInteger(iColumn)) {
3807          assert(fabs(value - solution[iColumn]) <= 1.0e-5);
3808          double value2 = floor(value + 0.5);
3809          if (dj[iColumn] < -1.0e-6) {
3810            // negative dj
3811            //assert (value2==colUpper[iColumn]);
3812            if (saveUpper[iColumn] == colUpper[iColumn]) {
3813              nAtUbNatural++;
3814              sortValue = 0.0;
3815              double value = -dj[iColumn];
3816              if (value > gap)
3817                nFix++;
3818              else if (gap < value * bestRatio)
3819                bestRatio = gap / value;
3820              if (saveLower[iColumn] != colLower[iColumn]) {
3821                nNotNeeded++;
3822                sortValue = 1.0e20;
3823              }
3824            } else if (saveLower[iColumn] == colUpper[iColumn]) {
3825              nAtLbFixed++;
3826              sortValue = dj[iColumn];
3827            } else {
3828              nAtOther++;
3829              sortValue = 0.0;
3830              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
3831                nNotNeeded++;
3832                sortValue = 1.0e20;
3833              }
3834            }
3835          } else if (dj[iColumn] > 1.0e-6) {
3836            // positive dj
3837            //assert (value2==colLower[iColumn]);
3838            if (saveLower[iColumn] == colLower[iColumn]) {
3839              nAtLbNatural++;
3840              sortValue = 0.0;
3841              double value = dj[iColumn];
3842              if (value > gap)
3843                nFix++;
3844              else if (gap < value * bestRatio)
3845                bestRatio = gap / value;
3846              if (saveUpper[iColumn] != colUpper[iColumn]) {
3847                nNotNeeded++;
3848                sortValue = 1.0e20;
3849              }
3850            } else if (saveUpper[iColumn] == colLower[iColumn]) {
3851              nAtUbFixed++;
3852              sortValue = -dj[iColumn];
3853            } else {
3854              nAtOther++;
3855              sortValue = 0.0;
3856              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
3857                nNotNeeded++;
3858                sortValue = 1.0e20;
3859              }
3860            }
3861          } else {
3862            // zero dj
3863            if (value2 == saveUpper[iColumn]) {
3864              nAtUbNaturalZero++;
3865              sortValue = 0.0;
3866              if (saveLower[iColumn] != colLower[iColumn]) {
3867                nNotNeeded++;
3868                sortValue = 1.0e20;
3869              }
3870            } else if (value2 == saveLower[iColumn]) {
3871              nAtLbNaturalZero++;
3872              sortValue = 0.0;
3873            } else {
3874              nAtOtherNatural++;
3875              sortValue = 0.0;
3876              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
3877                nNotNeeded++;
3878                sortValue = 1.0e20;
3879              }
3880            }
3881          }
3882#if HOTSTART == 3
3883          sortValue = -fabs(dj[iColumn]);
3884#endif
3885        }
3886        hotstartSolution_[iColumn] = value;
3887        saveLower[iColumn] = sortValue;
3888        order[iColumn] = iColumn;
3889      }
3890      COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3891        nFix, bestRatio, gap));
3892      int nNeg = 0;
3893      CoinSort_2(saveLower, saveLower + numberColumns, order);
3894      for (int i = 0; i < numberColumns; i++) {
3895        if (saveLower[i] < 0.0) {
3896          nNeg++;
3897#if HOTSTART == 2 || HOTSTART == 3
3898          // swap sign ?
3899          saveLower[i] = -saveLower[i];
3900#endif
3901        }
3902      }
3903      CoinSort_2(saveLower, saveLower + nNeg, order);
3904      for (int i = 0; i < numberColumns; i++) {
3905#if HOTSTART == 1
3906        hotstartPriorities_[order[i]] = 100;
3907#else
3908          hotstartPriorities_[order[i]] = -(i + 1);
3909#endif
3910      }
3911      COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3912        nAtLbNatural,
3913        nAtUbNatural,
3914        nAtLbNaturalZero,
3915        nAtUbNaturalZero,
3916        nAtLbFixed,
3917        nAtUbFixed,
3918        nAtOther,
3919        nAtOtherNatural, nNotNeeded, nNeg));
3920      delete[] saveLower;
3921      delete[] saveUpper;
3922      if (!saveCompare) {
3923        // create depth first comparison
3924        saveCompare = nodeCompare_;
3925        // depth first
3926        nodeCompare_ = new CbcCompareDepth();
3927        tree_->setComparison(*nodeCompare_);
3928      }
3929    }
3930#endif
3931    newNode = new CbcNode;
3932    // Set objective value (not so obvious if NLP etc)
3933    setObjectiveValue(newNode, NULL);
3934    anyAction = -1;
3935    // To make depth available we may need a fake node
3936    CbcNode fakeNode;
3937    if (!currentNode_) {
3938      // Not true if sub trees assert (!numberNodes_);
3939      currentNode_ = &fakeNode;
3940    }
3941    phase_ = 3;
3942    // only allow 1000 passes
3943    int numberPassesLeft = 1000;
3944    // This is first crude step
3945    if (numberAnalyzeIterations_ && !parentModel_) {
3946      delete[] analyzeResults_;
3947      //int numberColumns = solver_->getNumCols();
3948      analyzeResults_ = new double[5 * numberIntegers_];
3949      numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3950      if (numberFixedAtRoot_ > 0) {
3951        COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3952        setPointers(solver_);
3953        numberFixedNow_ = numberFixedAtRoot_;
3954      } else if (numberFixedAtRoot_ < 0) {
3955        COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3956        anyAction = -2;
3957        delete newNode;
3958        newNode = NULL;
3959        feasible = false;
3960      }
3961    }
3962    OsiSolverBranch *branches = NULL;
3963    if (feasible)
3964      anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3965        NULL, NULL, NULL, branches);
3966    if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3967      if (anyAction != -2) {
3968        // zap parent nodeInfo
3969#ifdef COIN_DEVELOP
3970        printf("zapping CbcNodeInfo %x\n", newNode->nodeInfo()->parent());
3971#endif
3972        if (newNode->nodeInfo())
3973          newNode->nodeInfo()->nullParent();
3974      }
3975      delete newNode;
3976      newNode = NULL;
3977      feasible = false;
3978    }
3979  }
3980  if (newNode && probingInfo_) {
3981    int number01 = probingInfo_->numberIntegers();
3982    //const fixEntry * entry = probingInfo_->fixEntries();
3983    const int *toZero = probingInfo_->toZero();
3984    //const int * toOne = probingInfo_->toOne();
3985    //const int * integerVariable = probingInfo_->integerVariable();
3986    if (toZero[number01]) {
3987      CglTreeProbingInfo info(*probingInfo_);
3988      if ((moreSpecialOptions2_ & 64) != 0 && !parentModel_) {
3989        /*
3990                Marginal idea. Further exploration probably good. Build some extra
3991                cliques from probing info. Not quite worth the effort?
3992              */
3993        CglProbing generator1;
3994        generator1.setUsingObjective(false);
3995        generator1.setMaxPass(1);
3996        generator1.setMaxPassRoot(1);
3997        generator1.setMaxLook(100);
3998        generator1.setRowCuts(3);
3999        generator1.setMaxElements(300);
4000        generator1.setMaxProbeRoot(solver_->getNumCols());
4001        CoinThreadRandom randomGenerator;
4002        //CglTreeProbingInfo info(solver_);
4003        info.level = 0;
4004        info.formulation_rows = solver_->getNumRows();
4005        info.inTree = false;
4006        if (parentModel_) {
4007          info.parentSolver = parentModel_->continuousSolver();
4008          // indicate if doing full search
4009          info.hasParent = ((specialOptions_ & 67108864) == 0) ? 1 : 2;
4010        } else {
4011          info.hasParent = 0;
4012          info.parentSolver = NULL;
4013        }
4014        info.originalColumns = originalColumns();
4015        info.randomNumberGenerator = &randomGenerator;
4016        info.pass = 4;
4017        generator1.setMode(8);
4018        OsiCuts cs;
4019        generator1.generateCutsAndModify(*solver_, cs, &info);
4020        // very clunky
4021        OsiSolverInterface *temp = generator1.cliqueModel(solver_, 2);
4022        CglPreProcess dummy;
4023        OsiSolverInterface *newSolver = dummy.cliqueIt(*temp, 0.0001);
4024        delete temp;
4025        OsiSolverInterface *fake = NULL;
4026        if (newSolver) {
4027#if 0
4028                int numberCliques = generator1.numberCliques();
4029                cliqueEntry * entry = generator1.cliqueEntry();
4030                cliqueType * type = new cliqueType [numberCliques];
4031                int * start = new int [numberCliques+1];
4032                start[numberCliques]=2*numberCliques;
4033                int n=0;
4034                for (int i=0;i<numberCliques;i++) {
4035                  start[i]=2*i;
4036                  setOneFixesInCliqueEntry(entry[2*i],true);
4037                  setOneFixesInCliqueEntry(entry[2*i+1],true);
4038                  type[i]=0;
4039                }
4040                fake = info.analyze(*solver_, 1,numberCliques,start,
4041                                    entry,type);
4042                delete [] type;
4043                delete [] entry;
4044#else
4045        fake = info.analyze(*newSolver, 1, -1);
4046#endif
4047          delete newSolver;
4048        } else {
4049          fake = info.analyze(*solver_, 1);
4050        }
4051        if (fake) {
4052          //fake->writeMps("fake");
4053          CglFakeClique cliqueGen(fake);
4054          cliqueGen.setStarCliqueReport(false);
4055          cliqueGen.setRowCliqueReport(false);
4056          cliqueGen.setMinViolation(0.1);
4057          addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
4058          generator_[numberCutGenerators_ - 1]->setTiming(true);
4059          for (int i = 0; i < numberCutGenerators_; i++) {
4060            CglKnapsackCover *cutGen = dynamic_cast< CglKnapsackCover * >(generator_[i]->generator());
4061            if (cutGen) {
4062              cutGen->createCliques(*fake, 2, 200, false);
4063            }
4064          }
4065        }
4066      }
4067      if (probingInfo_->packDown()) {
4068#if CBC_USEFUL_PRINTING > 1
4069        printf("%d implications on %d 0-1\n", toZero[number01], number01);
4070#endif
4071        // Create a cut generator that remembers implications discovered at root.
4072        CglImplication implication(probingInfo_);
4073        addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
4074        generator_[numberCutGenerators_ - 1]->setGlobalCuts(true);
4075        generator_[numberCutGenerators_ - 1]->setTiming(true);
4076      } else {
4077        delete probingInfo_;
4078        probingInfo_ = NULL;
4079      }
4080    } else {
4081      delete probingInfo_;
4082
4083      probingInfo_ = NULL;
4084    }
4085  }
4086  /*
4087      At this point, the root subproblem is infeasible or fathomed by bound
4088      (newNode == NULL), or we're live with an objective value that satisfies the
4089      current objective cutoff.
4090    */
4091  assert(!newNode || newNode->objectiveValue() <= cutoff);
4092  // Save address of root node as we don't want to delete it
4093  /*
4094      The common case is that the lp relaxation is feasible but doesn't satisfy
4095      integrality (i.e., newNode->branchingObject(), indicating we've been able to
4096      select a branching variable). Remove any cuts that have gone slack due to
4097      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
4098      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
4099      addCuts()). If, by some miracle, we have an integral solution at the root
4100      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
4101      a valid solution for use by setBestSolution().
4102    */
4103  CoinWarmStartBasis *lastws = NULL;
4104  if (feasible && newNode->branchingObject()) {
4105    if (resolved) {
4106      takeOffCuts(cuts, false, NULL);
4107#ifdef CHECK_CUT_COUNTS
4108      {
4109        printf("Number of rows after chooseBranch fix (root)"
4110               "(active only) %d\n",
4111          numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_);
4112        const CoinWarmStartBasis *debugws = dynamic_cast< const CoinWarmStartBasis * >(solver_->getWarmStart());
4113        debugws->print();
4114        delete debugws;
4115      }
4116#endif
4117    }
4118    //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
4119    //newNode->nodeInfo()->addCuts(cuts,
4120    //                   newNode->numberBranches(),whichGenerator_) ;
4121    if (lastws)
4122      delete lastws;
4123    lastws = dynamic_cast< CoinWarmStartBasis * >(solver_->getWarmStart());
4124  }
4125  /*
4126      Continuous data to be used later
4127    */
4128  continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4129  continuousInfeasibilities_ = 0;
4130  if (newNode) {
4131    continuousObjective_ = newNode->objectiveValue();
4132    delete[] continuousSolution_;
4133    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4134      numberColumns);
4135    continuousInfeasibilities_ = newNode->numberUnsatisfied();
4136  }
4137  /*
4138      Bound may have changed so reset in objects
4139    */
4140  {
4141    int i;
4142    for (i = 0; i < numberObjects_; i++)
4143      object_[i]->resetBounds(solver_);
4144  }
4145  /*
4146      Feasible? Then we should have either a live node prepped for future
4147      expansion (indicated by variable() >= 0), or (miracle of miracles) an
4148      integral solution at the root node.
4149
4150      initializeInfo sets the reference counts in the nodeInfo object.  Since
4151      this node is still live, push it onto the heap that holds the live set.
4152    */
4153  if (newNode) {
4154    if (newNode->branchingObject()) {
4155      newNode->initializeInfo();
4156      tree_->push(newNode);
4157      // save pointer to root node - so can pick up bounds
4158      if (!topOfTree_)
4159        topOfTree_ = dynamic_cast< CbcFullNodeInfo * >(newNode->nodeInfo());
4160      if (statistics_) {
4161        if (numberNodes2_ == maximumStatistics_) {
4162          maximumStatistics_ = 2 * maximumStatistics_;
4163          CbcStatistics **temp = new CbcStatistics *[maximumStatistics_];
4164          memset(temp, 0, maximumStatistics_ * sizeof(CbcStatistics *));
4165          memcpy(temp, statistics_, numberNodes2_ * sizeof(CbcStatistics *));
4166          delete[] statistics_;
4167          statistics_ = temp;
4168        }
4169        assert(!statistics_[numberNodes2_]);
4170        statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4171      }
4172      numberNodes2_++;
4173#ifdef CHECK_NODE
4174      printf("Node %x on tree\n", newNode);
4175#endif
4176    } else {
4177      // continuous is integer
4178      double objectiveValue = newNode->objectiveValue();
4179      setBestSolution(CBC_SOLUTION, objectiveValue,
4180        solver_->getColSolution());
4181      if (eventHandler) {
4182        // we are stopping anyway so no need to test return code
4183        eventHandler->event(CbcEventHandler::solution);
4184      }
4185      delete newNode;
4186      newNode = NULL;
4187    }
4188  }
4189
4190  if (printFrequency_ <= 0) {
4191    printFrequency_ = 1000;
4192    if (getNumCols() > 2000)
4193      printFrequency_ = 100;
4194  }
4195  /*
4196      It is possible that strong branching fixes one variable and then the code goes round
4197      again and again.  This can take too long.  So we need to warn user - just once.
4198    */
4199  numberLongStrong_ = 0;
4200  CbcNode *createdNode = NULL;
4201#ifdef CBC_THREAD
4202  if ((specialOptions_ & 2048) != 0)
4203    numberThreads_ = 0;
4204  if (numberThreads_) {
4205    nodeCompare_->sayThreaded(); // need to use addresses
4206    master_ = new CbcBaseModel(*this,
4207      (parallelMode() < -1) ? 1 : 0);
4208    masterThread_ = master_->masterThread();
4209  }
4210#endif
4211#ifdef COIN_HAS_CLP
4212  {
4213    OsiClpSolverInterface *clpSolver
4214      = dynamic_cast< OsiClpSolverInterface * >(solver_);
4215    if (clpSolver && !parentModel_) {
4216      clpSolver->computeLargestAway();
4217    }
4218  }
4219#endif
4220  /*
4221      At last, the actual branch-and-cut search loop, which will iterate until
4222      the live set is empty or we hit some limit (integrality gap, time, node
4223      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4224      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4225      add the node to the live set.
4226
4227      The first action is to winnow the live set to remove nodes which are worse
4228      than the current objective cutoff.
4229    */
4230  if (solver_->getRowCutDebuggerAlways()) {
4231    OsiRowCutDebugger *debuggerX = const_cast< OsiRowCutDebugger * >(solver_->getRowCutDebuggerAlways());
4232    const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger();
4233    if (!debugger) {
4234      // infeasible!!
4235      printf("before search\n");
4236      const double *lower = solver_->getColLower();
4237      const double *upper = solver_->getColUpper();
4238      const double *solution = debuggerX->optimalSolution();
4239      int numberColumns = solver_->getNumCols();
4240      for (int i = 0; i < numberColumns; i++) {
4241        if (solver_->isInteger(i)) {
4242          if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
4243            printf("**** ");
4244          printf("%d %g <= %g <= %g\n",
4245            i, lower[i], solution[i], upper[i]);
4246        }
4247      }
4248      //abort();
4249    }
4250  }
4251  {
4252    // may be able to change cutoff now
4253    double cutoff = getCutoff();
4254    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
4255    if (cutoff > bestObjective_ - increment) {
4256      cutoff = bestObjective_ - increment;
4257      setCutoff(cutoff);
4258    }
4259  }
4260#ifdef CBC_THREAD
4261  bool goneParallel = false;
4262#endif
4263#define MAX_DEL_NODE 1
4264  CbcNode *delNode[MAX_DEL_NODE + 1];
4265  int nDeleteNode = 0;
4266  // For Printing etc when parallel
4267  int lastEvery1000 = 0;
4268  int lastPrintEvery = 0;
4269  int numberConsecutiveInfeasible = 0;
4270#define PERTURB_IN_FATHOM
4271#ifdef PERTURB_IN_FATHOM
4272  // allow in fathom
4273  if ((moreSpecialOptions_ & 262144) != 0)
4274    specialOptions_ |= 131072;
4275#endif
4276  while (true) {
4277    lockThread();
4278#ifdef COIN_HAS_CLP
4279    // See if we want dantzig row choice
4280    goToDantzig(100, savePivotMethod);
4281#endif
4282    //#define REPORT_DYNAMIC 2
4283#if REPORT_DYNAMIC
4284    if (numberNodes_ && !parentModel_ && (tree_->empty() || (numberNodes_ % 10000) == 0)) {
4285      // Get average up and down costs
4286      double averageUp = 0.0;
4287      double averageDown = 0.0;
4288      int numberUp = 0;
4289      int numberDown = 0;
4290      int minTimesDown = COIN_INT_MAX;
4291      int maxTimesDown = 0;
4292      int neverBranchedDown = 0;
4293      int infeasibleTimesDown = 0;
4294      int minTimesUp = COIN_INT_MAX;
4295      int maxTimesUp = 0;
4296      int infeasibleTimesUp = 0;
4297      int neverBranchedUp = 0;
4298      int neverBranched = 0;
4299      int i;
4300      int numberInts = 0;
4301      bool endOfSearch = tree_->empty();
4302      int numberUp2 = 0;
4303      int numberDown2 = 0;
4304      for (i = 0; i < numberObjects_; i++) {
4305        OsiObject *object = object_[i];
4306        CbcSimpleIntegerDynamicPseudoCost *dynamicObject = dynamic_cast< CbcSimpleIntegerDynamicPseudoCost * >(object);
4307        if (dynamicObject) {
4308          numberInts++;
4309          if (dynamicObject->numberTimesUp() || dynamicObject->numberTimesDown()) {
4310            int nUp = 0;
4311            int nDown = 0;
4312            double up = 0.0;
4313            double down = 0.0;
4314            if (dynamicObject->numberTimesUp()) {
4315              numberUp++;
4316              nUp = dynamicObject->numberTimesUp();
4317              minTimesUp = CoinMin(minTimesUp, nUp);
4318              maxTimesUp = CoinMax(maxTimesUp, nUp);
4319              up = dynamicObject->upDynamicPseudoCost();
4320              averageUp += up;
4321              numberUp2 += nUp;
4322              infeasibleTimesUp += dynamicObject->numberTimesUpInfeasible();
4323            } else {
4324              neverBranchedUp++;
4325            }
4326            if (dynamicObject->numberTimesDown()) {
4327              numberDown++;
4328              nDown = dynamicObject->numberTimesDown();
4329              minTimesDown = CoinMin(minTimesDown, nDown);
4330              maxTimesDown = CoinMax(maxTimesDown, nDown);
4331              down = dynamicObject->downDynamicPseudoCost();
4332              averageDown += down;
4333              numberDown2 += dynamicObject->numberTimesDown();
4334              infeasibleTimesDown += dynamicObject->numberTimesDownInfeasible();
4335            } else {
4336              neverBranchedDown++;
4337            }
4338#if REPORT_DYNAMIC > 1
4339#if REPORT_DYNAMIC == 2
4340            if (endOfSearch && numberIntegers_ < 400) {
4341#elif REPORT_DYNAMIC == 3
4342            if (endOfSearch) {
4343#else
4344            {
4345#endif
4346              dynamicObject->print(0, 0.0);
4347            }
4348#endif
4349          } else {
4350            neverBranched++;
4351#if REPORT_DYNAMIC > 2
4352#if REPORT_DYNAMIC == 3
4353            if (endOfSearch && numberIntegers_ < 400) {
4354#elif REPORT_DYNAMIC == 4
4355            if (endOfSearch) {
4356#else
4357            {
4358#endif
4359              printf("col %d - never branched on\n", dynamicObject->columnNumber());
4360            }
4361#endif
4362          }
4363        }
4364      }
4365      if (numberUp)
4366        averageUp /= static_cast< double >(numberUp);
4367      else
4368        averageUp = 0.0;
4369      if (numberDown)
4370        averageDown /= static_cast< double >(numberDown);
4371      else
4372        averageDown = 0.0;
4373      printf("Report for %d variables (%d never branched on) after %d nodes - total solves down %d up %d\n",
4374        numberInts, neverBranched, numberNodes_, numberDown2, numberUp2);
4375      if ((neverBranchedDown || neverBranchedUp) && endOfSearch)
4376        printf("odd %d never branched down and %d never branched up\n",
4377          neverBranchedDown, neverBranchedUp);
4378      printf("down average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
4379        static_cast< double >(numberDown2) / numberDown, infeasibleTimesDown, averageDown,
4380        minTimesDown, maxTimesDown);
4381      printf("up average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
4382        static_cast< double >(numberUp2) / numberUp, infeasibleTimesUp, averageUp,
4383        minTimesUp, maxTimesUp);
4384    }
4385#endif
4386    if (tree_->empty()) {
4387#ifdef CBC_THREAD
4388      if (parallelMode() > 0 && master_) {
4389        int anyLeft = master_->waitForThreadsInTree(0);
4390        if (!anyLeft) {
4391          master_->stopThreads(-1);
4392          break;
4393        }
4394      } else {
4395        break;
4396      }
4397#else
4398    break;
4399#endif
4400    } else {
4401      unlockThread();
4402    }
4403    // If done 50/100 nodes see if worth trying reduction
4404    if (numberNodes_ >= nextCheckRestart) {
4405      if (nextCheckRestart < 100)
4406        nextCheckRestart = 100;
4407      else
4408        nextCheckRestart = COIN_INT_MAX;
4409#ifdef COIN_HAS_CLP
4410      OsiClpSolverInterface *clpSolver
4411        = dynamic_cast< OsiClpSolverInterface * >(solver_);
4412      if (clpSolver && ((specialOptions_ & 131072) == 0) && true) {
4413        ClpSimplex *simplex = clpSolver->getModelPtr();
4414        int perturbation = simplex->perturbation();
4415#if CBC_USEFUL_PRINTING > 1
4416        printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4417          numberIterations_, saveNumberIterations,
4418          numberSolves_, saveNumberSolves, perturbation);
4419#endif
4420        if (perturbation == 50 && (numberIterations_ - saveNumberIterations) < 8 * (numberSolves_ - saveNumberSolves)) {
4421          // switch off perturbation
4422          simplex->setPerturbation(100);
4423#if CBC_USEFUL_PRINTING > 1
4424          printf("Perturbation switched off\n");
4425#endif
4426        }
4427      }
4428#endif
4429      /*
4430              Decide if we want to do a restart.
4431            */
4432      if (saveSolver && (specialOptions_ & (512 + 32768)) != 0) {
4433        bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() && (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4434        int numberColumns = getNumCols();
4435        if (tryNewSearch) {
4436          // adding increment back allows current best - tiny bit weaker
4437          checkCutoffForRestart = getCutoff() + getCutoffIncrement();
4438#if CBC_USEFUL_PRINTING > 1
4439          printf("after %d nodes, cutoff %g - looking\n",
4440            numberNodes_, getCutoff());
4441#endif
4442          saveSolver->resolve();
4443          double direction = saveSolver->getObjSense();
4444          double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction;
4445          double tolerance;
4446          saveSolver->getDblParam(OsiDualTolerance, tolerance);
4447          if (gap <= 0.0)
4448            gap = tolerance;
4449          gap += 100.0 * tolerance;
4450          double integerTolerance = getDblParam(CbcIntegerTolerance);
4451
4452          const double *lower = saveSolver->getColLower();
4453          const double *upper = saveSolver->getColUpper();
4454          const double *solution = saveSolver->getColSolution();
4455          const double *reducedCost = saveSolver->getReducedCost();
4456
4457          int numberFixed = 0;
4458          int numberFixed2 = 0;
4459#ifdef COIN_DEVELOP
4460          printf("gap %g\n", gap);
4461#endif
4462          for (int i = 0; i < numberIntegers_; i++) {
4463            int iColumn = integerVariable_[i];
4464            double djValue = direction * reducedCost[iColumn];
4465            if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4466              if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4467                //printf("%d to lb on dj of %g - bounds %g %g\n",
4468                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4469                saveSolver->setColUpper(iColumn, lower[iColumn]);
4470                numberFixed++;
4471              } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4472                //printf("%d to ub on dj of %g - bounds %g %g\n",
4473                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4474                saveSolver->setColLower(iColumn, upper[iColumn]);
4475                numberFixed++;
4476              }
4477            } else {
4478              //printf("%d has dj of %g - already fixed to %g\n",
4479              //     iColumn,djValue,lower[iColumn]);
4480              numberFixed2++;
4481            }
4482          }
4483#ifdef COIN_DEVELOP
4484          if ((specialOptions_ & 1) != 0) {
4485            const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger();
4486            if (debugger) {
4487              printf("Contains optimal\n");
4488              OsiSolverInterface *temp = saveSolver->clone();
4489              const double *solution = debugger->optimalSolution();
4490              const double *lower = temp->getColLower();
4491              const double *upper = temp->getColUpper();
4492              int n = temp->getNumCols();
4493              for (int i = 0; i < n; i++) {
4494                if (temp->isInteger(i)) {
4495                  double value = floor(solution[i] + 0.5);
4496                  assert(value >= lower[i] && value <= upper[i]);
4497                  temp->setColLower(i, value);
4498                  temp->setColUpper(i, value);
4499                }
4500              }
4501              temp->writeMps("reduced_fix");
4502              delete temp;
4503              saveSolver->writeMps("reduced");
4504            } else {
4505              abort();
4506            }
4507          }
4508          printf("Restart could fix %d integers (%d already fixed)\n",
4509            numberFixed + numberFixed2, numberFixed2);
4510#endif
4511          numberFixed += numberFixed2;
4512          if (numberFixed * 10 < numberColumns && numberFixed * 4 < numberIntegers_)
4513            tryNewSearch = false;
4514        }
4515#ifdef CONFLICT_CUTS
4516        // temporary
4517        //if ((moreSpecialOptions_&4194304)!=0)
4518        //tryNewSearch=false;
4519#endif
4520        if (tryNewSearch) {
4521          // back to solver without cuts?
4522          OsiSolverInterface *solver2 = saveSolver->clone();
4523          const double *lower = saveSolver->getColLower();
4524          const double *upper = saveSolver->getColUpper();
4525          for (int i = 0; i < numberIntegers_; i++) {
4526            int iColumn = integerVariable_[i];
4527            solver2->setColLower(iColumn, lower[iColumn]);
4528            solver2->setColUpper(iColumn, upper[iColumn]);
4529          }
4530          // swap
4531          delete saveSolver;
4532          saveSolver = solver2;
4533          double *newSolution = new double[numberColumns];
4534          double objectiveValue = checkCutoffForRestart;
4535          // Save the best solution so far.
4536          CbcSerendipity heuristic(*this);
4537          if (bestSolution_)
4538            heuristic.setInputSolution(bestSolution_, bestObjective_);
4539          // Magic number
4540          heuristic.setFractionSmall(0.8);
4541          // `pumpTune' to stand-alone solver for explanations.
4542          heuristic.setFeasibilityPumpOptions(1008013);
4543          // Use numberNodes to say how many are original rows
4544          heuristic.setNumberNodes(continuousSolver_->getNumRows());
4545#ifdef COIN_DEVELOP
4546          if (continuousSolver_->getNumRows() < solver_->getNumRows())
4547            printf("%d rows added ZZZZZ\n",
4548              solver_->getNumRows() - continuousSolver_->getNumRows());
4549#endif
4550          int returnCode = heuristic.smallBranchAndBound(saveSolver,
4551            -1, newSolution,
4552            objectiveValue,
4553            checkCutoffForRestart, "Reduce");
4554          if (returnCode < 0) {
4555#ifdef COIN_DEVELOP
4556            printf("Restart - not small enough to do search after fixing\n");
4557#endif
4558            delete[] newSolution;
4559          } else {
4560            // 1 for sol'n, 2 for finished, 3 for both
4561            if ((returnCode & 1) != 0) {
4562              // increment number of solutions so other heuristics can test
4563              numberSolutions_++;
4564              numberHeuristicSolutions_++;
4565              lastHeuristic_ = NULL;
4566              setBestSolution(CBC_ROUNDING, objectiveValue, newSolution);
4567            }
4568            delete[] newSolution;
4569#ifdef CBC_THREAD
4570            if (master_) {
4571              lockThread();
4572              if (parallelMode() > 0) {
4573                while (master_->waitForThreadsInTree(0)) {
4574                  lockThread();
4575                  double dummyBest;
4576                  tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4577                  //unlockThread();
4578                }
4579              } else {
4580                double dummyBest;
4581                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4582              }
4583              master_->waitForThreadsInTree(2);
4584              delete master_;
4585              master_ = NULL;
4586              masterThread_ = NULL;
4587            }
4588#endif
4589            if (tree_->size()) {
4590              double dummyBest;
4591              tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4592            }
4593            break;
4594          }
4595        }
4596        delete saveSolver;
4597        saveSolver = NULL;
4598      }
4599    }
4600    /*
4601          Check for abort on limits: node count, solution count, time, integrality gap.
4602        */
4603    if (!(numberNodes_ < intParam_[CbcMaxNumNode] && numberSolutions_ < intParam_[CbcMaxNumSol] && !maximumSecondsReached() && !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))) {
4604      // out of loop
4605      break;
4606    }
4607#ifdef BONMIN
4608    assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4609#endif
4610// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4611// we need to replace the comparison function to dive, and that requires reordering to retain the
4612// heap property.
4613#define DIVE_WHEN 1000
4614#define DIVE_STOP 2000
4615    int kNode = numberNodes_ % 4000;
4616    if (numberNodes_ < 100000 && kNode > DIVE_WHEN && kNode <= DIVE_STOP) {
4617      if (!parallelMode()) {
4618        if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4619          CbcCompareDefault *compare = dynamic_cast< CbcCompareDefault * >(nodeCompare_);
4620          // Don't interfere if user has replaced the compare function.
4621          if (compare) {
4622            //printf("Redoing tree\n");
4623            compare->startDive(this);
4624            numberConsecutiveInfeasible = 0;
4625          }
4626        }
4627      }
4628    }
4629    // replace current cutoff?
4630    if (cutoff > getCutoff()) {
4631      double newCutoff = getCutoff();
4632      if (analyzeResults_) {
4633        // see if we could fix any (more)
4634        int n = 0;
4635        double *newLower = analyzeResults_;
4636        double *objLower = newLower + numberIntegers_;
4637        double *newUpper = objLower + numberIntegers_;
4638        double *objUpper = newUpper + numberIntegers_;
4639        for (int i = 0; i < numberIntegers_; i++) {
4640          if (objLower[i] > newCutoff) {
4641            n++;
4642            if (objUpper[i] > newCutoff) {
4643              newCutoff = -COIN_DBL_MAX;
4644              break;
4645            }
4646            // add as global cut
4647            objLower[i] = -COIN_DBL_MAX;
4648            OsiRowCut rc;
4649            rc.setLb(newLower[i]);
4650            rc.setUb(COIN_DBL_MAX);
4651            double one = 1.0;
4652            rc.setRow(1, integerVariable_ + i, &one, false);
4653            rc.setGloballyValidAsInteger(2);
4654            globalCuts_.addCutIfNotDuplicate(rc);
4655          } else if (objUpper[i] > newCutoff) {
4656            n++;
4657            // add as global cut
4658            objUpper[i] = -COIN_DBL_MAX;
4659            OsiRowCut rc;
4660            rc.setLb(-COIN_DBL_MAX);
4661            rc.setUb(newUpper[i]);
4662            double one = 1.0;
4663            rc.setRow(1, integerVariable_ + i, &one, false);
4664            rc.setGloballyValidAsInteger(2);
4665            globalCuts_.addCutIfNotDuplicate(rc);
4666          }
4667        }
4668        if (newCutoff == -COIN_DBL_MAX) {
4669          COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4670        } else if (n > numberFixedNow_) {
4671          COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4672          numberFixedNow_ = n;
4673        }
4674      }
4675      if (eventHandler) {
4676        if (!eventHandler->event(CbcEventHandler::solution)) {
4677          eventHappened_ = true; // exit
4678        }
4679        newCutoff = getCutoff();
4680      }
4681      lockThread();
4682      /*
4683              Clean the tree to reflect the new solution, then see if the
4684              node comparison predicate wants to make any changes. If so,
4685              call setComparison for the side effect of rebuilding the heap.
4686            */
4687      tree_->cleanTree(this, newCutoff, bestPossibleObjective_);
4688      if (nodeCompare_->newSolution(this) || nodeCompare_->newSolution(this, continuousObjective_, continuousInfeasibilities_)) {
4689        tree_->setComparison(*nodeCompare_);
4690      }
4691      if (tree_->empty()) {
4692        continue;
4693      }
4694      unlockThread();
4695    }
4696    cutoff = getCutoff();
4697    /*
4698            Periodic activities: Opportunities to
4699            + tweak the nodeCompare criteria,
4700            + check if we've closed the integrality gap enough to quit,
4701            + print a summary line to let the user know we're working
4702        */
4703    if (numberNodes_ >= lastEvery1000) {
4704      lockThread();
4705#ifdef COIN_HAS_CLP
4706      // See if we want dantzig row choice
4707      goToDantzig(1000, savePivotMethod);
4708#endif
4709      lastEvery1000 = numberNodes_ + 1000;
4710      bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_);
4711#ifdef CHECK_CUT_SIZE
4712      verifyCutSize(tree_, *this);
4713#endif
4714      // redo tree if requested
4715      if (redoTree)
4716        tree_->setComparison(*nodeCompare_);
4717      unlockThread();
4718    }
4719    // Had hotstart before, now switched off
4720    if (saveCompare && !hotstartSolution_) {
4721      // hotstart switched off
4722      delete nodeCompare_; // off depth first
4723      nodeCompare_ = saveCompare;
4724      saveCompare = NULL;
4725      // redo tree
4726      lockThread();
4727      tree_->setComparison(*nodeCompare_);
4728      unlockThread();
4729    }
4730    if (numberNodes_ >= lastPrintEvery) {
4731      lastPrintEvery = numberNodes_ + printFrequency_;
4732      lockThread();
4733      int nNodes = tree_->size();
4734
4735      //MODIF PIERRE
4736      bestPossibleObjective_ = tree_->getBestPossibleObjective();
4737#ifdef CBC_THREAD
4738      if (parallelMode() > 0 && master_) {
4739        // need to adjust for ones not on tree
4740        int numberThreads = master_->numberThreads();
4741        for (int i = 0; i < numberThreads; i++) {
4742          CbcThread *child = master_->child(i);
4743          if (child->node()) {
4744            // adjust
4745            double value = child->node()->objectiveValue();
4746            bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
4747          }
4748        }
4749      }
4750#endif
4751      unlockThread();
4752#if CBC_USEFUL_PRINTING > 1
4753      if (getCutoff() < 1.0e20) {
4754        if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 && !parentModel_)
4755          printf("model cutoff in status %g, best %g, increment %g\n",
4756            getCutoff(), bestObjective_, getCutoffIncrement());
4757        assert(getCutoff() < bestObjective_ - getCutoffIncrement() + 1.0e-6 + 1.0e-10 * fabs(bestObjective_));
4758      }
4759#endif
4760      if (!intParam_[CbcPrinting]) {
4761        // Parallel may not have any nodes
4762        if (!nNodes)
4763          bestPossibleObjective_ = lastBestPossibleObjective;
4764        else
4765          lastBestPossibleObjective = bestPossibleObjective_;
4766        messageHandler()->message(CBC_STATUS, messages())
4767          << numberNodes_ << CoinMax(nNodes, 1) << bestObjective_ << bestPossibleObjective_
4768          << getCurrentSeconds()
4769          << CoinMessageEol;
4770      } else if (intParam_[CbcPrinting] == 1) {
4771        messageHandler()->message(CBC_STATUS2, messages())
4772          << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4773          << tree_->lastDepth() << tree_->lastUnsatisfied()
4774          << tree_->lastObjective() << numberIterations_
4775          << getCurrentSeconds()
4776          << CoinMessageEol;
4777      } else if (!numberExtraIterations_) {
4778        messageHandler()->message(CBC_STATUS2, messages())
4779          << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4780          << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4781          << getCurrentSeconds()
4782          << CoinMessageEol;
4783      } else {
4784        messageHandler()->message(CBC_STATUS3, messages())
4785          << numberNodes_ << numberFathoms_ << numberExtraNodes_ << nNodes
4786          << bestObjective_ << bestPossibleObjective_
4787          << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4788          << getCurrentSeconds()
4789          << CoinMessageEol;
4790      }
4791#ifdef COIN_HAS_NTY
4792      if (symmetryInfo_)
4793        symmetryInfo_->statsOrbits(this, 1);
4794#endif
4795#if PRINT_CONFLICT == 1
4796      if (numberConflictCuts > lastNumberConflictCuts) {
4797        double length = lengthConflictCuts / numberConflictCuts;
4798        printf("%d new conflict cuts - total %d - average length %g\n",
4799          numberConflictCuts - lastNumberConflictCuts,
4800          numberConflictCuts, length);
4801        lastNumberConflictCuts = numberConflictCuts;
4802      }
4803#endif
4804      if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4805        eventHappened_ = true; // exit
4806      }
4807    }
4808    // See if can stop on gap
4809    if (canStopOnGap()) {
4810      stoppedOnGap_ = true;
4811    }
4812
4813#ifdef CHECK_NODE_FULL
4814    verifyTreeNodes(tree_, *this);
4815#endif
4816#ifdef CHECK_CUT_COUNTS
4817    verifyCutCounts(tree_, *this);
4818#endif
4819    /*
4820          Now we come to the meat of the loop. To create the active subproblem, we'll
4821          pop the most promising node in the live set, rebuild the subproblem it
4822          represents, and then execute the current arm of the branch to create the
4823          active subproblem.
4824        */
4825    CbcNode *node = NULL;
4826#ifdef CBC_THREAD
4827    if (!parallelMode() || parallelMode() == -1) {
4828#endif
4829      node = tree_->bestNode(cutoff);
4830      // Possible one on tree worse than cutoff
4831      // Weird comparison function can leave ineligible nodes on tree
4832      if (!node || node->objectiveValue() > cutoff)
4833        continue;
4834      // Do main work of solving node here
4835      doOneNode(this, node, createdNode);
4836#ifdef JJF_ZERO
4837      if (node) {
4838        if (createdNode) {
4839          printf("Node %d depth %d, created %d depth %d\n",
4840            node->nodeNumber(), node->depth(),
4841            createdNode->nodeNumber(), createdNode->depth());
4842        } else {
4843          printf("Node %d depth %d,  no created node\n",
4844            node->nodeNumber(), node->depth());
4845        }
4846      } else if (createdNode) {
4847        printf("Node exhausted, created %d depth %d\n",
4848          createdNode->nodeNumber(), createdNode->depth());
4849      } else {
4850        printf("Node exhausted,  no created node\n");
4851        numberConsecutiveInfeasible = 2;
4852      }
4853#endif
4854      //if (createdNode)
4855      //numberConsecutiveInfeasible=0;
4856      //else
4857      //numberConsecutiveInfeasible++;
4858#ifdef CBC_THREAD
4859    } else if (parallelMode() > 0) {
4860      //lockThread();
4861      //node = tree_->bestNode(cutoff) ;
4862      // Possible one on tree worse than cutoff
4863      if (true || !node || node->objectiveValue() > cutoff) {
4864        assert(master_);
4865        if (master_) {
4866          int anyLeft = master_->waitForThreadsInTree(1);
4867          // may need to go round again
4868          if (anyLeft) {
4869            continue;
4870          } else {
4871            master_->stopThreads(-1);
4872          }
4873        }
4874      }
4875      //unlockThread();
4876    } else {
4877      // Deterministic parallel
4878      if ((tree_->size() < CoinMax(numberThreads_, 8) || hotstartSolution_) && !goneParallel) {
4879        node = tree_->bestNode(cutoff);
4880        // Possible one on tree worse than cutoff
4881        if (!node || node->objectiveValue() > cutoff)
4882          continue;
4883        // Do main work of solving node here
4884        doOneNode(this, node, createdNode);
4885        assert(createdNode);
4886        if (!createdNode->active()) {
4887          delete createdNode;
4888          createdNode = NULL;
4889        } else {
4890          // Say one more pointing to this
4891          node->nodeInfo()->increment();
4892          tree_->push(createdNode);
4893        }
4894        if (node->active()) {
4895          assert(node->nodeInfo());
4896          if (node->nodeInfo()->numberBranchesLeft()) {
4897            tree_->push(node);
4898          } else {
4899            node->setActive(false);
4900          }
4901        } else {
4902          if (node->nodeInfo()) {
4903            if (!node->nodeInfo()->numberBranchesLeft())
4904              node->nodeInfo()->allBranchesGone(); // can clean up
4905            // So will delete underlying stuff
4906            node->setActive(true);
4907          }
4908          delNode[nDeleteNode++] = node;
4909          node = NULL;
4910        }
4911        if (nDeleteNode >= MAX_DEL_NODE) {
4912          for (int i = 0; i < nDeleteNode; i++) {
4913            //printf("trying to del %d %x\n",i,delNode[i]);
4914            delete delNode[i];
4915            //printf("done to del %d %x\n",i,delNode[i]);
4916          }
4917          nDeleteNode = 0;
4918        }
4919      } else {
4920        // Split and solve
4921        master_->deterministicParallel();
4922        goneParallel = true;
4923      }
4924    }
4925#endif
4926  }
4927  if (nDeleteNode) {
4928    for (int i = 0; i < nDeleteNode; i++) {
4929      delete delNode[i];
4930    }
4931    nDeleteNode = 0;
4932  }
4933#ifdef CBC_THREAD
4934  if (master_) {
4935    master_->stopThreads(-1);
4936    master_->waitForThreadsInTree(2);
4937    // adjust time to allow for children on some systems
4938    //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4939  }
4940#endif
4941  /*
4942      End of the non-abort actions. The next block of code is executed if we've
4943      aborted because we hit one of the limits. Clean up by deleting the live set
4944      and break out of the node processing loop. Note that on an abort, node may
4945      have been pushed back onto the tree for further processing, in which case
4946      it'll be deleted in cleanTree. We need to check.
4947    */
4948  if (!(numberNodes_ < intParam_[CbcMaxNumNode] && numberSolutions_ < intParam_[CbcMaxNumSol] && !maximumSecondsReached() && !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))) {
4949    if (tree_->size()) {
4950      double dummyBest;
4951      tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
4952#if 0 // Does not seem to be needed def CBC_THREAD
4953            if (parallelMode() > 0 && master_) {
4954              // see if any dangling nodes
4955              int numberThreads = master_->numberThreads();
4956              for (int i=0;i<numberThreads;i++) {
4957                CbcThread * child = master_->child(i);
4958                //if (child->createdNode())
4959                //printf("CHILD_NODE %p\n",child->createdNode());
4960                delete child->createdNode();
4961              }
4962            }
4963#endif
4964    }
4965    delete nextRowCut_;
4966    /* order is important here:
4967         * maximumSecondsReached() should be checked before eventHappened_ and
4968         * isNodeLimitReached() should be checked after eventHappened_
4969         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4970         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4971         */
4972    if (stoppedOnGap_) {
4973      messageHandler()->message(CBC_GAP, messages())
4974        << bestObjective_ - bestPossibleObjective_
4975        << dblParam_[CbcAllowableGap]
4976        << dblParam_[CbcAllowableFractionGap] * 100.0
4977        << CoinMessageEol;
4978      secondaryStatus_ = 2;
4979      status_ = 0;
4980    } else if (maximumSecondsReached()) {
4981      handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol;
4982      secondaryStatus_ = 4;
4983      status_ = 1;
4984    } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4985      handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol;
4986      secondaryStatus_ = 6;
4987      status_ = 1;
4988    } else if (isNodeLimitReached()) {
4989      handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol;
4990      secondaryStatus_ = 3;
4991      status_ = 1;
4992    } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4993      handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol;
4994      secondaryStatus_ = 8;
4995      status_ = 1;
4996    } else {
4997      handler_->message(CBC_EVENT, messages_) << CoinMessageEol;
4998      secondaryStatus_ = 5;
4999      status_ = 5;
5000    }
5001  }
5002#ifdef CBC_THREAD
5003  if (master_) {
5004    delete master_;
5005    master_ = NULL;
5006    masterThread_ = NULL;
5007  }
5008#endif
5009  /*
5010      That's it, we've exhausted the search tree, or broken out of the loop because
5011      we hit some limit on evaluation.
5012
5013      We may have got an intelligent tree so give it one more chance
5014    */
5015  // Tell solver we are not in Branch and Cut
5016  solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL);
5017  tree_->endSearch();
5018  //  If we did any sub trees - did we give up on any?
5019  if (numberStoppedSubTrees_)
5020    status_ = 1;
5021  numberNodes_ += numberExtraNodes_;
5022  numberIterations_ += numberExtraIterations_;
5023  if (eventHandler) {
5024    eventHandler->event(CbcEventHandler::endSearch);
5025  }
5026  if (!status_) {
5027    // Set best possible unless stopped on gap
5028    if (secondaryStatus_ != 2)
5029      bestPossibleObjective_ = bestObjective_;
5030    handler_->message(CBC_END_GOOD, messages_)
5031      << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
5032      << CoinMessageEol;
5033  } else {
5034    handler_->message(CBC_END, messages_)
5035      << bestObjective_ << bestPossibleObjective_
5036      << numberIterations_ << numberNodes_ << getCurrentSeconds()
5037      << CoinMessageEol;
5038  }
5039  if ((moreSpecialOptions_ & 4194304) != 0) {
5040    // Conflict cuts
5041    int numberCuts = globalCuts_.sizeRowCuts();
5042    int nConflict = 0;
5043    double sizeConflict = 0.0;
5044    for (int i = 0; i < numberCuts; i++) {
5045      OsiRowCut2 *cut = globalCuts_.cut(i);
5046      if (cut->whichRow() == 1) {
5047        nConflict++;
5048        sizeConflict += cut->row().getNumElements();
5049      }
5050    }
5051    if (nConflict) {
5052      sizeConflict /= nConflict;
5053      char general[200];
5054      sprintf(general, "%d conflict cuts generated - average length %g",
5055        nConflict, sizeConflict);
5056      messageHandler()->message(CBC_GENERAL,
5057        messages())
5058        << general << CoinMessageEol;
5059    }
5060  }
5061  if (numberStrongIterations_)
5062    handler_->message(CBC_STRONG_STATS, messages_)
5063      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
5064      << strongInfo_[1] << CoinMessageEol;
5065  if (!numberExtraNodes_)
5066    handler_->message(CBC_OTHER_STATS, messages_)
5067      << maximumDepthActual_
5068      << numberDJFixed_ << CoinMessageEol;
5069  else
5070    handler_->message(CBC_OTHER_STATS2, messages_)
5071      << maximumDepthActual_
5072      << numberDJFixed_ << numberFathoms_ << numberExtraNodes_ << numberExtraIterations_
5073      << CoinMessageEol;
5074#ifdef COIN_HAS_NTY
5075  if (symmetryInfo_)
5076    symmetryInfo_->statsOrbits(this, 1);
5077#endif
5078  if (doStatistics == 100) {
5079    for (int i = 0; i < numberObjects_; i++) {
5080      CbcSimpleIntegerDynamicPseudoCost *obj = dynamic_cast< CbcSimpleIntegerDynamicPseudoCost * >(object_[i]);
5081      if (obj)
5082        obj->print();
5083    }
5084  }
5085  if (statistics_) {
5086    // report in some way
5087    int *lookup = new int[numberObjects_];
5088    int i;
5089    for (i = 0; i < numberObjects_; i++)
5090      lookup[i] = -1;
5091    bool goodIds = false; //true;
5092    for (i = 0; i < numberObjects_; i++) {
5093      int iColumn = object_[i]->columnNumber();
5094      if (iColumn >= 0 && iColumn < numberColumns) {
5095        if (lookup[i] == -1) {
5096          lookup[i] = iColumn;
5097        } else {
5098          goodIds = false;
5099          break;
5100        }
5101      } else {
5102        goodIds = false;
5103        break;
5104      }
5105    }
5106    if (!goodIds) {
5107      delete[] lookup;
5108      lookup = NULL;
5109    }
5110    if (doStatistics >= 3) {
5111      printf("  node parent depth column   value                    obj      inf\n");
5112      for (i = 0; i < numberNodes2_; i++) {
5113        statistics_[i]->print(lookup);
5114      }
5115    }
5116    if (doStatistics > 1) {
5117      // Find last solution
5118      int k;
5119      for (k = numberNodes2_ - 1; k >= 0; k--) {
5120        if (statistics_[k]->endingObjective() != COIN_DBL_MAX && !statistics_[k]->endingInfeasibility())
5121          break;
5122      }
5123      if (k >= 0) {
5124        int depth = statistics_[k]->depth();
5125        int *which = new int[depth + 1];
5126        for (i = depth; i >= 0; i--) {
5127          which[i] = k;
5128          k = statistics_[k]->parentNode();
5129        }
5130        printf("  node parent depth column   value                    obj      inf\n");
5131        for (i = 0; i <= depth; i++) {
5132          statistics_[which[i]]->print(lookup);
5133        }
5134        delete[] which;
5135      }
5136    }
5137    // now summary
5138    int maxDepth = 0;
5139    double averageSolutionDepth = 0.0;
5140    int numberSolutions = 0;
5141    double averageCutoffDepth = 0.0;
5142    double averageSolvedDepth = 0.0;
5143    int numberCutoff = 0;
5144    int numberDown = 0;
5145    int numberFirstDown = 0;
5146    double averageInfDown = 0.0;
5147    double averageObjDown = 0.0;
5148    int numberCutoffDown = 0;
5149    int numberUp = 0;
5150    int numberFirstUp = 0;
5151    double averageInfUp = 0.0;
5152    double averageObjUp = 0.0;
5153    int numberCutoffUp = 0;
5154    double averageNumberIterations1 = 0.0;
5155    double averageValue = 0.0;
5156    for (i = 0; i < numberNodes2_; i++) {
5157      int depth = statistics_[i]->depth();
5158      int way = statistics_[i]->way();
5159      double value = statistics_[i]->value();
5160      double startingObjective = statistics_[i]->startingObjective();
5161      int startingInfeasibility = statistics_[i]->startingInfeasibility();
5162      double endingObjective = statistics_[i]->endingObjective();
5163      int endingInfeasibility = statistics_[i]->endingInfeasibility();
5164      maxDepth = CoinMax(depth, maxDepth);
5165      // Only for completed
5166      averageNumberIterations1 += statistics_[i]->numberIterations();
5167      averageValue += value;
5168      if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
5169        numberSolutions++;
5170        averageSolutionDepth += depth;
5171      }
5172      if (endingObjective == COIN_DBL_MAX) {
5173        numberCutoff++;
5174        averageCutoffDepth += depth;
5175        if (way < 0) {
5176          numberDown++;
5177          numberCutoffDown++;
5178          if (way == -1)
5179            numberFirstDown++;
5180        } else {
5181          numberUp++;
5182          numberCutoffUp++;
5183          if (way == 1)
5184            numberFirstUp++;
5185        }
5186      } else {
5187        averageSolvedDepth += depth;
5188        if (way < 0) {
5189          numberDown++;
5190          averageInfDown += startingInfeasibility - endingInfeasibility;
5191          averageObjDown += endingObjective - startingObjective;
5192          if (way == -1)
5193            numberFirstDown++;
5194        } else {
5195          numberUp++;
5196          averageInfUp += startingInfeasibility - endingInfeasibility;
5197          averageObjUp += endingObjective - startingObjective;
5198          if (way == 1)
5199            numberFirstUp++;
5200        }
5201      }
5202    }
5203    // Now print
5204    if (numberSolutions)
5205      averageSolutionDepth /= static_cast< double >(numberSolutions);
5206    int numberSolved = numberNodes2_ - numberCutoff;
5207    double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
5208      - numberIterationsAtContinuous;
5209    if (numberCutoff) {
5210      averageCutoffDepth /= static_cast< double >(numberCutoff);
5211      averageNumberIterations2 /= static_cast< double >(numberCutoff);
5212    }
5213    if (numberNodes2_)
5214      averageValue /= static_cast< double >(numberNodes2_);
5215    if (numberSolved) {
5216      averageNumberIterations1 /= static_cast< double >(numberSolved);
5217      averageSolvedDepth /= static_cast< double >(numberSolved);
5218    }
5219    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
5220      numberSolutions, averageSolutionDepth);
5221    printf("average value of variable being branched on was %g\n",
5222      averageValue);
5223    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
5224      numberCutoff, averageCutoffDepth, averageNumberIterations2);
5225    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
5226      numberSolved, averageSolvedDepth, averageNumberIterations1);
5227    if (numberDown) {
5228      averageInfDown /= static_cast< double >(numberDown);
5229      averageObjDown /= static_cast< double >(numberDown);
5230    }
5231    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5232      numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
5233      averageInfDown, averageObjDown);
5234    if (numberUp) {
5235      averageInfUp /= static_cast< double >(numberUp);
5236      averageObjUp /= static_cast< double >(numberUp);
5237    }
5238    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5239      numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
5240      averageInfUp, averageObjUp);
5241    for (i = 0; i < numberNodes2_; i++)
5242      delete statistics_[i];
5243    delete[] statistics_;
5244    statistics_ = NULL;
5245    maximumStatistics_ = 0;
5246    delete[] lookup;
5247  }
5248  /*
5249      If we think we have a solution, restore and confirm it with a call to
5250      setBestSolution().  We need to reset the cutoff value so as not to fathom
5251      the solution on bounds.  Note that calling setBestSolution( ..., true)
5252      leaves the continuousSolver_ bounds vectors fixed at the solution value.
5253
5254      Running resolve() here is a failsafe --- setBestSolution has already
5255      reoptimised using the continuousSolver_. If for some reason we fail to
5256      prove optimality, run the problem again after instructing the solver to
5257      tell us more.
5258
5259      If all looks good, replace solver_ with continuousSolver_, so that the
5260      outside world will be able to obtain information about the solution using
5261      public methods.
5262
5263      Don't replace if we are trying to save cuts
5264    */
5265  if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && ((specialOptions_ & 8388608) == 0 || (specialOptions_ & 2048) != 0)) {
5266    setCutoff(1.0e50); // As best solution should be worse than cutoff
5267    // change cutoff as constraint if wanted
5268    if (cutoffRowNumber_ >= 0) {
5269      if (solver_->getNumRows() > cutoffRowNumber_)
5270        solver_->setRowUpper(cutoffRowNumber_, 1.0e50);
5271    }
5272    // also in continuousSolver_
5273    if (continuousSolver_) {
5274      // Solvers know about direction
5275      double direction = solver_->getObjSense();
5276      continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
5277    }
5278    phase_ = 5;
5279    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
5280    if ((specialOptions_ & 4) == 0)
5281      bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
5282    setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1);
5283    continuousSolver_->resolve();
5284    // Deal with funny variables
5285    if ((moreSpecialOptions2_ & 32768) != 0)
5286      cleanBounds(continuousSolver_, NULL);
5287    if (!continuousSolver_->isProvenOptimal()) {
5288      continuousSolver_->messageHandler()->setLogLevel(2);
5289      continuousSolver_->initialSolve();
5290    }
5291    delete solver_;
5292    // above deletes solverCharacteristics_
5293    solverCharacteristics_ = NULL;
5294    solver_ = continuousSolver_;
5295    setPointers(solver_);
5296    continuousSolver_ = NULL;
5297  }
5298  /*
5299      Clean up dangling objects. continuousSolver_ may already be toast.
5300    */
5301  delete lastws;
5302  if (saveObjects) {
5303    for (int i = 0; i < numberObjects_; i++)
5304      delete saveObjects[i];
5305    delete[] saveObjects;
5306  }
5307  numberStrong_ = saveNumberStrong;
5308  numberBeforeTrust_ = saveNumberBeforeTrust;
5309  delete[] whichGenerator_;
5310  whichGenerator_ = NULL;
5311  delete[] lowerBefore;
5312  delete[] upperBefore;
5313  delete[] walkback_;
5314  walkback_ = NULL;
5315  delete[] lastNodeInfo_;
5316  lastNodeInfo_ = NULL;
5317  delete[] lastNumberCuts_;
5318  lastNumberCuts_ = NULL;
5319  delete[] lastCut_;
5320  lastCut_ = NULL;
5321  delete[] addedCuts_;
5322  addedCuts_ = NULL;
5323  //delete persistentInfo;
5324  // Get rid of characteristics
5325  solverCharacteristics_ = NULL;
5326  if (continuousSolver_) {
5327    delete continuousSolver_;
5328    continuousSolver_ = NULL;
5329  }
5330  /*
5331      Destroy global cuts by replacing with an empty OsiCuts object.
5332    */
5333  globalCuts_ = CbcRowCuts();
5334  delete globalConflictCuts_;
5335  globalConflictCuts_ = NULL;
5336  if (!bestSolution_ && (specialOptions_ & 8388608) == 0 && false) {
5337    // make sure lp solver is infeasible
5338    int numberColumns = solver_->getNumCols();
5339    const double *columnLower = solver_->getColLower();
5340    int iColumn;
5341    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5342      if (solver_->isInteger(iColumn))
5343        solver_->setColUpper(iColumn, columnLower[iColumn]);
5344    }
5345    solver_->initialSolve();
5346  }
5347#ifdef COIN_HAS_CLP
5348  {
5349    OsiClpSolverInterface *clpSolver
5350      = dynamic_cast< OsiClpSolverInterface * >(solver_);
5351    if (clpSolver) {
5352      // Possible restore of pivot method
5353      if (savePivotMethod) {
5354        // model may have changed
5355        savePivotMethod->setModel(NULL);
5356        clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
5357        delete savePivotMethod;
5358      }
5359      clpSolver->setLargestAway(-1.0);
5360    }
5361  }
5362#endif
5363  if ((fastNodeDepth_ >= 1000000 || (moreSpecialOptions_ & 33554432) != 0)
5364    && !parentModel_) {
5365    // delete object off end
5366    delete object_[numberObjects_];
5367    if ((moreSpecialOptions_ & 33554432) == 0)
5368      fastNodeDepth_ -= 1000000;
5369  }
5370  delete saveSolver;
5371  // Undo preprocessing performed during BaB.
5372  if (strategy_ && strategy_->preProcessState() > 0) {
5373    // undo preprocessing
5374    CglPreProcess *process = strategy_->process();
5375    assert(process);
5376    int n = originalSolver->getNumCols();
5377    if (bestSolution_) {
5378      delete[] bestSolution_;
5379      bestSolution_ = new double[n];
5380      process->postProcess(*solver_);
5381    }
5382    strategy_->deletePreProcess();
5383    // Solution now back in originalSolver
5384    delete solver_;
5385    solver_ = originalSolver;
5386    if (bestSolution_) {
5387      bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
5388      memcpy(bestSolution_, solver_->getColSolution(), n * sizeof(double));
5389    }
5390    // put back original objects if there were any
5391    if (originalObject) {
5392      int iColumn;
5393      assert(ownObjects_);
5394      for (iColumn = 0; iColumn < numberObjects_; iColumn++)
5395        delete object_[iColumn];
5396      delete[] object_;
5397      numberObjects_ = numberOriginalObjects;
5398      object_ = originalObject;
5399      delete[] integerVariable_;
5400      numberIntegers_ = 0;
5401      for (iColumn = 0; iColumn < n; iColumn++) {
5402        if (solver_->isInteger(iColumn))
5403          numberIntegers_++;
5404      }
5405      integerVariable_ = new int[numberIntegers_];
5406      numberIntegers_ = 0;
5407      for (iColumn = 0; iColumn < n; iColumn++) {
5408        if (solver_->isInteger(iColumn))
5409          integerVariable_[numberIntegers_++] = iColumn;
5410      }
5411    }
5412  }
5413  if (flipObjective)
5414    flipModel();
5415#ifdef COIN_HAS_CLP
5416  {
5417    OsiClpSolverInterface *clpSolver
5418      = dynamic_cast< OsiClpSolverInterface * >(solver_);
5419    if (clpSolver)
5420      clpSolver->setFakeObjective(reinterpret_cast< double * >(NULL));
5421  }
5422#endif
5423  moreSpecialOptions_ = saveMoreSpecialOptions;
5424  return;
5425}
5426
5427// Solve the initial LP relaxation
5428void CbcModel::initialSolve()
5429{
5430  assert(solver_);
5431  // Double check optimization directions line up
5432  dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
5433  // Check if bounds are all integral (as may get messed up later)
5434  checkModel();
5435  if (!solverCharacteristics_) {
5436    OsiBabSolver *solverCharacteristics = dynamic_cast< OsiBabSolver * >(solver_->getAuxiliaryInfo());
5437    if (solverCharacteristics) {
5438      solverCharacteristics_ = solverCharacteristics;
5439    } else {
5440      // replace in solver
5441      OsiBabSolver defaultC;
5442      solver_->setAuxiliaryInfo(&defaultC);
5443      solverCharacteristics_ = dynamic_cast< OsiBabSolver * >(solver_->getAuxiliaryInfo());
5444    }
5445  }
5446  solverCharacteristics_->setSolver