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

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

don't allow restart if off objects

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