source: trunk/Cbc/src/CbcModel.cpp

Last change on this file was 2698, checked in by stefan, 2 weeks ago

do memcpy only if source is not NULL

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