source: branches/devel/Cbc/src/CbcModel.cpp @ 642

Last change on this file since 642 was 642, checked in by forrest, 12 years ago

update branches/devel for threads

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 349.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include "CbcConfig.h"
9
10#include <string>
11//#define CBC_DEBUG 1
12//#define CHECK_CUT_COUNTS
13//#define CHECK_NODE_FULL
14//#define NODE_LOG
15//#define GLOBAL_CUTS_JUST_POINTERS
16#include <cassert>
17#include <cmath>
18#include <cfloat>
19
20#ifdef COIN_HAS_CLP
21// include Presolve from Clp
22#include "ClpPresolve.hpp"
23#include "OsiClpSolverInterface.hpp"
24#endif
25
26#include "CbcEventHandler.hpp"
27
28#include "OsiSolverInterface.hpp"
29#include "OsiAuxInfo.hpp"
30#include "OsiSolverBranch.hpp"
31#include "OsiChooseVariable.hpp"
32#include "CoinWarmStartBasis.hpp"
33#include "CoinPackedMatrix.hpp"
34#include "CoinHelperFunctions.hpp"
35#include "CbcBranchActual.hpp"
36#include "CbcBranchDynamic.hpp"
37#include "CbcHeuristic.hpp"
38#include "CbcHeuristicFPump.hpp"
39#include "CbcModel.hpp"
40#include "CbcTreeLocal.hpp"
41#include "CbcStatistics.hpp"
42#include "CbcStrategy.hpp"
43#include "CbcMessage.hpp"
44#include "OsiRowCut.hpp"
45#include "OsiColCut.hpp"
46#include "OsiRowCutDebugger.hpp"
47#include "OsiCuts.hpp"
48#include "CbcCountRowCut.hpp"
49#include "CbcCutGenerator.hpp"
50#include "CbcFeasibilityBase.hpp"
51#include "CbcFathom.hpp"
52// include Probing
53#include "CglProbing.hpp"
54// include preprocessing
55#include "CglPreProcess.hpp"
56#include "CglDuplicateRow.hpp"
57#include "CglStored.hpp"
58
59#include "CoinTime.hpp"
60#include "CoinMpsIO.hpp"
61
62#include "CbcCompareActual.hpp"
63#include "CbcTree.hpp"
64//#define CBC_THREAD
65#ifdef CBC_THREAD
66#include <pthread.h>
67//#define CBC_THREAD_DEBUG 1
68// To Pass across to doOneNode
69typedef struct {
70  CbcModel * baseModel;
71  CbcModel * thisModel;
72  CbcNode * node; // filled in every time
73  CbcNode * createdNode; // filled in every time on return
74  pthread_t threadIdOfBase;
75  pthread_mutex_t * mutex; // for locking data
76  pthread_mutex_t * mutex2; // for waking up threads
77  pthread_cond_t * condition2; // for waking up thread
78  int returnCode; // -1 available, 0 busy, 1 finished , 2??
79  double timeLocked;
80  double timeWaitingToLock;
81  double timeWaitingToStart;
82  double timeInThread;
83  int numberTimesLocked;
84  int numberTimesUnlocked;
85  int numberTimesWaitingToStart;
86  int saveStuff[2];
87  struct timespec absTime;
88  bool locked;
89#if CBC_THREAD_DEBUG
90  int threadNumber;
91#endif
92} threadStruct;
93static void * doNodesThread(void * voidInfo);
94static void * doCutsThread(void * voidInfo);
95#endif
96/* Various functions local to CbcModel.cpp */
97
98namespace {
99
100//-------------------------------------------------------------------
101// Returns the greatest common denominator of two
102// positive integers, a and b, found using Euclid's algorithm
103//-------------------------------------------------------------------
104static int gcd(int a, int b) 
105{
106  int remainder = -1;
107  // make sure a<=b (will always remain so)
108  if(a > b) {
109    // Swap a and b
110    int temp = a;
111    a = b;
112    b = temp;
113  }
114  // if zero then gcd is nonzero (zero may occur in rhs of packed)
115  if (!a) {
116    if (b) {
117      return b;
118    } else {
119      printf("**** gcd given two zeros!!\n");
120      abort();
121    }
122  }
123  while (remainder) {
124    remainder = b % a;
125    b = a;
126    a = remainder;
127  }
128  return b;
129}
130
131
132
133#ifdef CHECK_NODE_FULL
134
135/*
136  Routine to verify that tree linkage is correct. The invariant that is tested
137  is
138
139  reference count = (number of actual references) + (number of branches left)
140
141  The routine builds a set of paired arrays, info and count, by traversing the
142  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
143  referenced (via the parent field) is recorded in count. Then a final check is
144  made to see if the numberPointingToThis_ field agrees.
145*/
146
147void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model)
148
149{if (model.getNodeCount()==661) return;  printf("*** CHECKING tree after %d nodes\n",model.getNodeCount()) ;
150 
151  int j ;
152  int nNodes = branchingTree->size() ;
153# define MAXINFO 1000
154  int *count = new int [MAXINFO] ;
155  CbcNodeInfo **info = new CbcNodeInfo*[MAXINFO] ;
156  int nInfo = 0 ;
157/*
158  Collect all CbcNodeInfo objects in info, by starting from each live node and
159  traversing back to the root. Nodes in the live set should have unexplored
160  branches remaining.
161
162  TODO: The `while (nodeInfo)' loop could be made to break on reaching a
163        common ancester (nodeInfo is found in info[k]). Alternatively, the
164        check could change to signal an error if nodeInfo is not found above a
165        common ancestor.
166*/
167  for (j = 0 ; j < nNodes ; j++)
168  { CbcNode *node = branchingTree->nodePointer(j) ;
169  if (!node)
170    continue;
171    CbcNodeInfo *nodeInfo = node->nodeInfo() ; 
172    int change = node->nodeInfo()->numberBranchesLeft() ;
173    assert(change) ;
174    while (nodeInfo)
175    { int k ;
176      for (k = 0 ; k < nInfo ; k++)
177      { if (nodeInfo == info[k]) break ; }
178      if (k == nInfo)
179      { assert(nInfo < MAXINFO) ;
180        nInfo++ ;
181        info[k] = nodeInfo ;
182        count[k] = 0 ; }
183      nodeInfo = nodeInfo->parent() ; } }
184/*
185  Walk the info array. For each nodeInfo, look up its parent in info and
186  increment the corresponding count.
187*/
188  for (j = 0 ; j < nInfo ; j++)
189  { CbcNodeInfo *nodeInfo = info[j] ;
190    nodeInfo = nodeInfo->parent() ;
191    if (nodeInfo)
192    { int k ;
193      for (k = 0 ; k < nInfo ; k++)
194      { if (nodeInfo == info[k]) break ; }
195      assert (k < nInfo) ;
196      count[k]++ ; } }
197/*
198  Walk the info array one more time and check that the invariant holds. The
199  number of references (numberPointingToThis()) should equal the sum of the
200  number of actual references (held in count[]) plus the number of potential
201  references (unexplored branches, numberBranchesLeft()).
202*/
203  for (j = 0;j < nInfo;j++) {
204    CbcNodeInfo * nodeInfo = info[j] ;
205    if (nodeInfo) {
206      int k ;
207      for (k = 0;k < nInfo;k++)
208        if (nodeInfo == info[k])
209          break ;
210      printf("Nodeinfo %x - %d left, %d count\n",
211             nodeInfo,
212             nodeInfo->numberBranchesLeft(),
213             nodeInfo->numberPointingToThis()) ;
214      assert(nodeInfo->numberPointingToThis() ==
215             count[k]+nodeInfo->numberBranchesLeft()) ; } }
216
217  delete [] count ;
218  delete [] info ;
219 
220  return ; }
221
222#endif  /* CHECK_NODE_FULL */
223
224
225
226#ifdef CHECK_CUT_COUNTS
227
228/*
229  Routine to verify that cut reference counts are correct.
230*/
231void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model)
232
233{ printf("*** CHECKING cuts after %d nodes\n",model.getNodeCount()) ;
234
235  int j ;
236  int nNodes = branchingTree->size() ;
237
238/*
239  cut.tempNumber_ exists for the purpose of doing this verification. Clear it
240  in all cuts. We traverse the tree by starting from each live node and working
241  back to the root. At each CbcNodeInfo, check for cuts.
242*/
243  for (j = 0 ; j < nNodes ; j++)
244  { CbcNode *node = branchingTree->nodePointer(j) ;
245    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
246    assert (node->nodeInfo()->numberBranchesLeft()) ;
247    while (nodeInfo)
248    { int k ;
249      for (k = 0 ; k < nodeInfo->numberCuts() ; k++)
250      { CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
251        if (cut) cut->tempNumber_ = 0; }
252      nodeInfo = nodeInfo->parent() ; } }
253/*
254  Walk the live set again, this time collecting the list of cuts in use at each
255  node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
256  that when we recreate the basis for a node, we compress out the slack cuts.
257*/
258  for (j = 0 ; j < nNodes ; j++)
259  { CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
260    CbcNode *node = branchingTree->nodePointer(j) ;
261    CbcNodeInfo *nodeInfo = node->nodeInfo(); 
262    int change = node->nodeInfo()->numberBranchesLeft() ;
263    printf("Node %d %x (info %x) var %d way %d obj %g",j,node,
264           node->nodeInfo(),node->columnNumber(),node->way(),
265           node->objectiveValue()) ;
266
267    model.addCuts1(node,debugws) ;
268
269    int i ;
270    int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
271    CbcCountRowCut **addedCuts = model.addedCuts() ;
272    for (i = 0 ; i < model.currentNumberCuts() ; i++)
273    { CoinWarmStartBasis::Status status = 
274        debugws->getArtifStatus(i+numberRowsAtContinuous) ;
275      if (status != CoinWarmStartBasis::basic && addedCuts[i])
276      { addedCuts[i]->tempNumber_ += change ; } }
277
278    while (nodeInfo)
279    { nodeInfo = nodeInfo->parent() ;
280      if (nodeInfo) printf(" -> %x",nodeInfo); }
281    printf("\n") ;
282    delete debugws ; }
283/*
284  The moment of truth: We've tallied up the references by direct scan of the
285  search tree. Check for agreement with the count in the cut.
286
287  TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
288*/
289  for (j = 0 ; j < nNodes ; j++)
290  { CbcNode *node = branchingTree->nodePointer(j) ;
291    CbcNodeInfo *nodeInfo = node->nodeInfo(); 
292    while (nodeInfo)
293    { int k ;
294      for (k = 0 ; k < nodeInfo->numberCuts() ; k++)
295      { CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
296        if (cut && cut->tempNumber_ >= 0)
297        { if (cut->tempNumber_ != cut->numberPointingToThis()) 
298            printf("mismatch %x %d %x %d %d\n",nodeInfo,k,
299                    cut,cut->tempNumber_,cut->numberPointingToThis()) ;
300          else
301            printf("   match %x %d %x %d %d\n", nodeInfo,k,
302                   cut,cut->tempNumber_,cut->numberPointingToThis()) ;
303          cut->tempNumber_ = -1 ; } }
304      nodeInfo = nodeInfo->parent() ; } }
305
306  return ; }
307
308#endif /* CHECK_CUT_COUNTS */
309
310
311//#define CHECK_CUT_SIZE
312#ifdef CHECK_CUT_SIZE
313
314/*
315  Routine to verify that cut reference counts are correct.
316*/
317void verifyCutSize (const CbcTree * branchingTree, CbcModel &model)
318{ 
319
320  int j ;
321  int nNodes = branchingTree->size() ;
322  int totalCuts=0;
323
324/*
325  cut.tempNumber_ exists for the purpose of doing this verification. Clear it
326  in all cuts. We traverse the tree by starting from each live node and working
327  back to the root. At each CbcNodeInfo, check for cuts.
328*/
329  for (j = 0 ; j < nNodes ; j++) {
330    CbcNode *node = branchingTree->nodePointer(j) ;
331    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
332    assert (node->nodeInfo()->numberBranchesLeft()) ;
333    while (nodeInfo) {
334      totalCuts += nodeInfo->numberCuts();
335      nodeInfo = nodeInfo->parent() ;
336    }
337  }
338  printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n",model.getNodeCount(),totalCuts) ;
339  return ;
340}
341
342#endif /* CHECK_CUT_SIZE */
343
344}
345
346 /* End unnamed namespace for CbcModel.cpp */
347
348
349static double trueIncrement=0.0;
350void 
351CbcModel::analyzeObjective ()
352/*
353  Try to find a minimum change in the objective function. The first scan
354  checks that there are no continuous variables with non-zero coefficients,
355  and grabs the largest objective coefficient associated with an unfixed
356  integer variable. The second scan attempts to scale up the objective
357  coefficients to a point where they are sufficiently close to integer that
358  we can pretend they are integer, and calculate a gcd over the coefficients
359  of interest. This will be the minimum increment for the scaled coefficients.
360  The final action is to scale the increment back for the original coefficients
361  and install it, if it's better than the existing value.
362
363  John's note: We could do better than this.
364
365  John's second note - apologies for changing s to z
366*/
367{ const double *objective = getObjCoefficients() ;
368  const double *lower = getColLower() ;
369  const double *upper = getColUpper() ;
370/*
371  Take a first scan to see if there are unfixed continuous variables in the
372  objective.  If so, the minimum objective change could be arbitrarily small.
373  Also pick off the maximum coefficient of an unfixed integer variable.
374
375  If the objective is found to contain only integer variables, set the
376  fathoming discipline to strict.
377*/
378  double maximumCost = 0.0 ;
379  trueIncrement=0.0;
380  bool possibleMultiple = true ;
381  int iColumn ;
382  int numberColumns = getNumCols() ;
383  for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
384  { if (upper[iColumn] > lower[iColumn]+1.0e-8)
385    { if (isInteger(iColumn)) 
386        maximumCost = CoinMax(maximumCost,fabs(objective[iColumn])) ;
387      else if (objective[iColumn]) 
388        possibleMultiple = false ; } }
389  setIntParam(CbcModel::CbcFathomDiscipline,possibleMultiple) ;
390/*
391  If a nontrivial increment is possible, try and figure it out. We're looking
392  for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
393  variables. Since the c<j> might not be integers, try and inflate them
394  sufficiently that they look like integers (and we'll deflate the gcd
395  later).
396
397  2520.0 is used as it is a nice multiple of 2,3,5,7
398*/
399    if (possibleMultiple&&maximumCost)
400    { int increment = 0 ;
401      double multiplier = 2520.0 ;
402      while (10.0*multiplier*maximumCost < 1.0e8)
403        multiplier *= 10.0 ;
404    int bigIntegers = 0; // Count of large costs which are integer
405    for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
406      if (upper[iColumn] > lower[iColumn]+1.0e-8) {
407        if (isInteger(iColumn)&&objective[iColumn]) {
408          double value = fabs(objective[iColumn])*multiplier ;
409          if (value <2.1e9) {
410            int nearest = (int) floor(value+0.5) ;
411            if (fabs(value-floor(value+0.5)) > 1.0e-8)
412              { increment = 0 ;
413              break ; }
414            else if (!increment)
415              { increment = nearest ; }
416            else
417              { increment = gcd(increment,nearest) ; }
418          } else {
419            // large value - may still be multiple of 1.0
420            value = fabs(objective[iColumn]);
421            if (fabs(value-floor(value+0.5)) > 1.0e-8) {
422              increment=0;
423              break;
424            } else {
425              bigIntegers++;
426            }
427          }
428        }
429      }
430    }
431
432/*
433  If the increment beats the current value for objective change, install it.
434*/
435      if (increment)
436      { double value = increment ;
437        double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ;
438        if (bigIntegers) {
439          // allow for 1.0
440          increment = gcd(increment,(int) multiplier);
441          value = increment;
442        }
443        value /= multiplier ;
444        trueIncrement=CoinMax(cutoff,value);;
445        if (value*0.999 > cutoff)
446        { messageHandler()->message(CBC_INTEGERINCREMENT,
447                                          messages())
448            << value << CoinMessageEol ;
449          setDblParam(CbcModel::CbcCutoffIncrement,value*0.999) ; } } }
450
451  return ; 
452}
453
454
455/**
456  \todo
457  Normally, it looks like we enter here from command dispatch in the main
458  routine, after calling the solver for an initial solution
459  (CbcModel::initialSolve, which simply calls the solver's initialSolve
460  routine.) The first thing we do is call resolve. Presumably there are
461  circumstances where this is nontrivial? There's also a call from
462  CbcModel::originalModel (tied up with integer presolve), which should be
463  checked.
464
465*/
466
467/*
468  The overall flow can be divided into three stages:
469    * Prep: Check that the lp relaxation remains feasible at the root. If so,
470      do all the setup for B&C.
471    * Process the root node: Generate cuts, apply heuristics, and in general do
472      the best we can to resolve the problem without B&C.
473    * Do B&C search until we hit a limit or exhaust the search tree.
474 
475  Keep in mind that in general there is no node in the search tree that
476  corresponds to the active subproblem. The active subproblem is represented
477  by the current state of the model,  of the solver, and of the constraint
478  system held by the solver.
479*/
480void CbcModel::branchAndBound(int doStatistics) 
481
482{
483/*
484  Capture a time stamp before we start.
485*/
486  dblParam_[CbcStartSeconds] = CoinCpuTime();
487  strongInfo_[0]=0;
488  strongInfo_[1]=0;
489  strongInfo_[2]=0;
490  numberStrongIterations_ = 0;
491#ifndef NDEBUG
492  {
493    int i;
494    int n = solver_->getNumCols();
495    const double *lower = solver_->getColLower() ;
496    const double *upper = solver_->getColUpper() ;
497    for (i=0;i<n;i++) {
498      assert (lower[i]<1.0e10);
499      assert (upper[i]>-1.0e10);
500    }
501    n = solver_->getNumRows();
502    lower = solver_->getRowLower() ;
503    upper = solver_->getRowUpper() ;
504    for (i=0;i<n;i++) {
505      assert (lower[i]<1.0e10);
506      assert (upper[i]>-1.0e10);
507    }
508  }
509#endif
510  // original solver (only set if pre-processing)
511  OsiSolverInterface * originalSolver=NULL;
512  int numberOriginalObjects=numberObjects_;
513  OsiObject ** originalObject = NULL;
514  // Set up strategies
515#if 0
516  std::string problemName ;
517  solver_->getStrParam(OsiProbName,problemName) ;
518  if (!strcmp(problemName.c_str(),"P0282")) solver_->activateRowCutDebugger(problemName.c_str()) ;
519#endif
520  if (strategy_) {
521    // May do preprocessing
522    originalSolver = solver_;
523    strategy_->setupOther(*this);
524    if (strategy_->preProcessState()) {
525      // pre-processing done
526      if (strategy_->preProcessState()<0) {
527        // infeasible
528        handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
529        status_ = 0 ;
530        secondaryStatus_ = 1;
531        originalContinuousObjective_ = COIN_DBL_MAX;
532        return ; 
533      } else if (numberObjects_&&object_) {
534        numberOriginalObjects=numberObjects_;
535        // redo sequence
536        numberIntegers_=0;
537        int numberColumns = getNumCols();
538        int nOrig = originalSolver->getNumCols();
539        CglPreProcess * process = strategy_->process();
540        assert (process);
541        const int * originalColumns = process->originalColumns();
542        // allow for cliques etc
543        nOrig = CoinMax(nOrig,originalColumns[numberColumns-1]+1);
544        // try and redo debugger
545        OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
546        if (debugger)
547          debugger->redoSolution(numberColumns,originalColumns);
548        originalObject = object_;
549        // object number or -1
550        int * temp = new int[nOrig];
551        int iColumn;
552        for (iColumn=0;iColumn<nOrig;iColumn++) 
553          temp[iColumn]=-1;
554        int iObject;
555        int nNonInt=0;
556        for (iObject=0;iObject<numberOriginalObjects;iObject++) {
557          iColumn = originalObject[iObject]->columnNumber();
558          if (iColumn<0) {
559            nNonInt++;
560          } else {
561            temp[iColumn]=iObject;
562          }
563        }
564        int numberNewIntegers=0;
565        int numberOldIntegers=0;
566        int numberOldOther=0;
567        for (iColumn=0;iColumn<numberColumns;iColumn++) {
568          int jColumn = originalColumns[iColumn];
569          if (temp[jColumn]>=0) {
570            int iObject= temp[jColumn];
571            CbcSimpleInteger * obj =
572              dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
573            if (obj) 
574              numberOldIntegers++;
575            else
576              numberOldOther++;
577          } else if (isInteger(iColumn)) {
578            numberNewIntegers++;
579          }
580        }
581        /*
582          Allocate an array to hold the indices of the integer variables.
583          Make a large enough array for all objects
584        */
585        numberObjects_= numberNewIntegers+numberOldIntegers+numberOldOther+nNonInt;
586        object_ = new OsiObject * [numberObjects_];
587        integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
588        /*
589          Walk the variables again, filling in the indices and creating objects for
590          the integer variables. Initially, the objects hold the index and upper &
591          lower bounds.
592        */
593        numberIntegers_=0;
594        for (iColumn=0;iColumn<numberColumns;iColumn++) {
595          int jColumn = originalColumns[iColumn];
596          if (temp[jColumn]>=0) {
597            int iObject= temp[jColumn];
598            CbcSimpleInteger * obj =
599              dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
600            if (obj) {
601              object_[numberIntegers_] = originalObject[iObject]->clone();
602              // redo ids etc
603              object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
604              integerVariable_[numberIntegers_++]=iColumn;
605            }
606          } else if (isInteger(iColumn)) {
607            object_[numberIntegers_] =
608              new CbcSimpleInteger(this,iColumn);
609            integerVariable_[numberIntegers_++]=iColumn;
610          }
611        }
612        numberObjects_=numberIntegers_;
613        // Now append other column stuff
614        for (iColumn=0;iColumn<numberColumns;iColumn++) {
615          int jColumn = originalColumns[iColumn];
616          if (temp[jColumn]>=0) {
617            int iObject= temp[jColumn];
618            CbcSimpleInteger * obj =
619              dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
620            if (!obj) {
621              object_[numberObjects_] = originalObject[iObject]->clone();
622              // redo ids etc
623              CbcObject * obj =
624              dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
625              assert (obj);
626              obj->redoSequenceEtc(this,numberColumns,originalColumns);
627              numberObjects_++;
628            }
629          }
630        }
631        // now append non column stuff
632        for (iObject=0;iObject<numberOriginalObjects;iObject++) {
633          iColumn = originalObject[iObject]->columnNumber();
634          if (iColumn<0) {
635            object_[numberObjects_] = originalObject[iObject]->clone();
636            // redo ids etc
637            CbcObject * obj =
638              dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
639            assert (obj);
640            obj->redoSequenceEtc(this,numberColumns,originalColumns);
641            numberObjects_++;
642          }
643        }
644        delete [] temp;
645        if (!numberObjects_)
646          handler_->message(CBC_NOINT,messages_) << CoinMessageEol ;
647      } else {
648        int numberColumns = getNumCols();
649        CglPreProcess * process = strategy_->process();
650        assert (process);
651        const int * originalColumns = process->originalColumns();
652        // try and redo debugger
653        OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
654        if (debugger)
655          debugger->redoSolution(numberColumns,originalColumns);
656      }
657    } else {
658      //no preprocessing
659      originalSolver=NULL;
660    }
661    strategy_->setupCutGenerators(*this);
662    strategy_->setupHeuristics(*this);
663    // Set strategy print level to models
664    strategy_->setupPrinting(*this,handler_->logLevel());
665  }
666  eventHappened_=false;
667  CbcEventHandler *eventHandler = getEventHandler() ;
668  // set up for probing
669  probingInfo_ = new CglTreeProbingInfo(solver_);
670
671  // Try for dominated columns
672  if ((specialOptions_&64)!=0) {
673    CglDuplicateRow dupcuts(solver_);
674    dupcuts.setMode(2);
675    CglStored * storedCuts = dupcuts.outDuplicates(solver_);
676    addCutGenerator(storedCuts,1,"StoredCuts from dominated");
677  }
678  if (!nodeCompare_)
679    nodeCompare_=new CbcCompareDefault();;
680  // See if hot start wanted
681  CbcCompareBase * saveCompare = NULL;
682  if (hotstartSolution_) {
683    if (strategy_&&strategy_->preProcessState()>0) {
684      CglPreProcess * process = strategy_->process();
685      assert (process);
686      int n = solver_->getNumCols();
687      const int * originalColumns = process->originalColumns();
688      // columns should be in order ... but
689      double * tempS = new double[n];
690      for (int i=0;i<n;i++) {
691        int iColumn = originalColumns[i];
692        tempS[i]=hotstartSolution_[iColumn];
693      }
694      delete [] hotstartSolution_;
695      hotstartSolution_=tempS;
696      if (hotstartPriorities_) {
697        int * tempP = new int [n];
698        for (int i=0;i<n;i++) {
699          int iColumn = originalColumns[i];
700          tempP[i]=hotstartPriorities_[iColumn];
701        }
702        delete [] hotstartPriorities_;
703        hotstartPriorities_=tempP;
704      }
705    }
706    saveCompare = nodeCompare_;
707    // depth first
708    nodeCompare_ = new CbcCompareDepth();
709  }
710  if (!problemFeasibility_)
711    problemFeasibility_=new CbcFeasibilityBase();
712# ifdef CBC_DEBUG
713  std::string problemName ;
714  solver_->getStrParam(OsiProbName,problemName) ;
715  printf("Problem name - %s\n",problemName.c_str()) ;
716  solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
717# endif
718/*
719  Assume we're done, and see if we're proven wrong.
720*/
721  status_ = 0 ;
722  secondaryStatus_ = 0;
723  phase_=0;
724/*
725  Scan the variables, noting the integer variables. Create an
726  CbcSimpleInteger object for each integer variable.
727*/
728  findIntegers(false) ;
729  // If dynamic pseudo costs then do
730  if (numberBeforeTrust_)
731    convertToDynamic();
732  // Set up char array to say if integer
733  delete [] integerInfo_;
734  {
735    int n = solver_->getNumCols();
736    integerInfo_ = new char [n];
737    for (int i=0;i<n;i++) {
738      if (solver_->isInteger(i))
739        integerInfo_[i]=1;
740      else
741        integerInfo_[i]=0;
742    }
743  }
744  if (preferredWay_) {
745    // set all unset ones
746    for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
747      CbcObject * obj =
748        dynamic_cast <CbcObject *>(object_[iObject]) ;
749      if (obj&&!obj->preferredWay())
750        obj->setPreferredWay(preferredWay_);
751    }
752  } 
753/*
754  Ensure that objects on the lists of OsiObjects, heuristics, and cut
755  generators attached to this model all refer to this model.
756*/
757  synchronizeModel() ;
758  if (!solverCharacteristics_) {
759    OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
760    if (solverCharacteristics) {
761      solverCharacteristics_ = solverCharacteristics;
762    } else {
763      // replace in solver
764      OsiBabSolver defaultC;
765      solver_->setAuxiliaryInfo(&defaultC);
766      solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
767    }
768  }
769  solverCharacteristics_->setSolver(solver_);
770  // Set so we can tell we are in initial phase in resolve
771  continuousObjective_ = -COIN_DBL_MAX ;
772/*
773  Solve the relaxation.
774
775  Apparently there are circumstances where this will be non-trivial --- i.e.,
776  we've done something since initialSolve that's trashed the solution to the
777  continuous relaxation.
778*/
779  bool feasible;
780  // If NLP then we assume already solved outside branchAndbound
781  if (!solverCharacteristics_->solverType()||solverCharacteristics_->solverType()==4) {
782    feasible=resolve(NULL,0) != 0 ;
783  } else {
784    // pick up given status
785    feasible = (solver_->isProvenOptimal() &&
786                !solver_->isDualObjectiveLimitReached()) ;
787  }
788  if (problemFeasibility_->feasible(this,0)<0) {
789    feasible=false; // pretend infeasible
790  }
791/*
792  If the linear relaxation of the root is infeasible, bail out now. Otherwise,
793  continue with processing the root node.
794*/
795  if (!feasible) {
796    status_ = 0 ;
797    if (!solver_->isProvenDualInfeasible()) {
798      handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
799      secondaryStatus_ = 1;
800    } else {
801      handler_->message(CBC_UNBOUNDED,messages_)<< CoinMessageEol ;
802      secondaryStatus_ = 7;
803    }
804    originalContinuousObjective_ = COIN_DBL_MAX;
805    solverCharacteristics_ = NULL;
806    return ;
807  }
808  // Convert to Osi if wanted
809  bool useOsiBranching=false;
810  //OsiBranchingInformation * persistentInfo = NULL;
811  if (branchingMethod_&&branchingMethod_->chooseMethod()) {
812    useOsiBranching=true;
813    //persistentInfo = new OsiBranchingInformation(solver_);
814    if (numberOriginalObjects) {
815      for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
816        CbcObject * obj =
817          dynamic_cast <CbcObject *>(object_[iObject]) ;
818        if (obj) {
819          CbcSimpleInteger * obj2 =
820            dynamic_cast <CbcSimpleInteger *>(obj) ;
821          if (obj2) {
822            // back to Osi land
823            object_[iObject]=obj2->osiObject();
824            delete obj;
825          } else {
826            OsiSimpleInteger * obj3 =
827              dynamic_cast <OsiSimpleInteger *>(obj) ;
828            if (!obj3) {
829              OsiSOS * obj4 =
830                dynamic_cast <OsiSOS *>(obj) ;
831              if (!obj4) {
832                CbcSOS * obj5 =
833                  dynamic_cast <CbcSOS *>(obj) ;
834                if (obj5) {
835                  // back to Osi land
836                  object_[iObject]=obj5->osiObject(solver_);
837                } else {
838                  printf("Code up CbcObject type in Osi land\n");
839                  abort();
840                }
841              }
842            }
843          }
844        }
845      }
846      // and add to solver if none
847      if (!solver_->numberObjects()) {
848        solver_->addObjects(numberObjects_,object_);
849      } else {
850        if (solver_->numberObjects()!=numberOriginalObjects) {
851          printf("should have trapped that solver has objects before\n");
852          abort();
853        }
854      }
855    } else {
856      // do from solver
857      deleteObjects(false);
858      solver_->findIntegersAndSOS(false);
859      numberObjects_=solver_->numberObjects();
860      object_ = solver_->objects();
861      ownObjects_ = false;
862    }
863    branchingMethod_->chooseMethod()->setSolver(solver_);
864  }
865  // take off heuristics if have to
866  if (numberHeuristics_) {
867    int numberOdd=0;
868    for (int i=0;i<numberObjects_;i++) {
869      if (!object_[i]->canDoHeuristics()) 
870        numberOdd++;
871    }
872    if (numberOdd) {
873      int k=0;
874      for (int i=0;i<numberHeuristics_;i++) {
875        if (!heuristic_[i]->canDealWithOdd())
876          delete heuristic_[i];
877        else
878          heuristic_[k++]=heuristic_[i];
879      }
880      if (!k) {
881        delete [] heuristic_;
882        heuristic_=NULL;
883      }
884      numberHeuristics_=k;
885      handler_->message(CBC_HEURISTICS_OFF,messages_)<< numberOdd<<CoinMessageEol ;
886    }
887  }
888  // Save objective (just so user can access it)
889  originalContinuousObjective_ = solver_->getObjValue();
890  bestPossibleObjective_=originalContinuousObjective_;
891  sumChangeObjective1_=0.0;
892  sumChangeObjective2_=0.0;
893/*
894  OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
895  Assuming it recognises the problem, when called upon it will check a cut to
896  see if it cuts off the optimal answer.
897*/
898  // If debugger exists set specialOptions_ bit
899  if (solver_->getRowCutDebuggerAlways())
900    specialOptions_ |= 1;
901
902# ifdef CBC_DEBUG
903  if ((specialOptions_&1)==0)
904    solver_->activateRowCutDebugger(problemName.c_str()) ;
905  if (solver_->getRowCutDebuggerAlways())
906    specialOptions_ |= 1;
907# endif
908
909/*
910  Begin setup to process a feasible root node.
911*/
912  bestObjective_ = CoinMin(bestObjective_,1.0e50) ;
913  numberSolutions_ = 0 ;
914  stateOfSearch_=0;
915  numberHeuristicSolutions_ = 0 ;
916  // Everything is minimization
917  { 
918    // needed to sync cutoffs
919    double value ;
920    solver_->getDblParam(OsiDualObjectiveLimit,value) ;
921    dblParam_[CbcCurrentCutoff]= value * solver_->getObjSense();
922  }
923  double cutoff=getCutoff() ;
924  double direction = solver_->getObjSense() ;
925  dblParam_[CbcOptimizationDirection]=direction;
926  if (cutoff < 1.0e20&&direction<0.0)
927    messageHandler()->message(CBC_CUTOFF_WARNING1,
928                                    messages())
929                                      << cutoff << -cutoff << CoinMessageEol ;
930  if (cutoff > bestObjective_)
931    cutoff = bestObjective_ ;
932  setCutoff(cutoff) ;
933/*
934  We probably already have a current solution, but just in case ...
935*/
936  int numberColumns = getNumCols() ;
937  if (!currentSolution_)
938    currentSolution_ = new double[numberColumns] ;
939  testSolution_ = currentSolution_;
940  /* Tell solver we are in Branch and Cut
941     Could use last parameter for subtle differences */
942  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
943#ifdef COIN_HAS_CLP
944 {
945   OsiClpSolverInterface * clpSolver
946     = dynamic_cast<OsiClpSolverInterface *> (solver_);
947   if (clpSolver) {
948     //#define CLP_QUICK_OPTIONS
949#ifdef CLP_QUICK_OPTIONS
950     // Try and re-use regions   
951     ClpSimplex * simplex = clpSolver->getModelPtr();
952     simplex->setPersistenceFlag(1);
953     clpSolver->deleteScaleFactors();
954     int value=131072;
955     clpSolver->setSpecialOptions(clpSolver->specialOptions()|value);
956     if ((clpSolver->specialOptions()&value)!=0)
957       simplex->setSpecialOptions(simplex->specialOptions()|value);
958     //if (simplex->numberRows()<50)
959     //simplex->setAlphaAccuracy(1.0);
960     //clpSolver->setSpecialOptions((clpSolver->specialOptions()&~128)|65536);
961#endif
962     if ((specialOptions_&32)==0) {
963       ClpSimplex * clpSimplex = clpSolver->getModelPtr();
964       // take off names
965       clpSimplex->dropNames();
966     }
967     // no crunch if mostly continuous
968     int numberColumns = solver_->getNumCols();
969     if (numberColumns>1000&&numberIntegers_*4<numberColumns)
970       clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~1));
971   }
972 }
973#endif
974/*
975  Create a copy of the solver, thus capturing the original (root node)
976  constraint system (aka the continuous system).
977*/
978  continuousSolver_ = solver_->clone() ;
979
980  numberRowsAtContinuous_ = getNumRows() ;
981/*
982  Check the objective to see if we can deduce a nontrivial increment. If
983  it's better than the current value for CbcCutoffIncrement, it'll be
984  installed.
985*/
986  if(solverCharacteristics_->reducedCostsAccurate())
987    analyzeObjective() ;
988/*
989  Set up for cut generation. addedCuts_ holds the cuts which are relevant for
990  the active subproblem. whichGenerator will be used to record the generator
991  that produced a given cut.
992*/
993  maximumWhich_ = 1000 ;
994  delete [] whichGenerator_;
995  whichGenerator_ = new int[maximumWhich_] ;
996  memset(whichGenerator_,0,maximumWhich_*sizeof(int));
997  maximumNumberCuts_ = 0 ;
998  currentNumberCuts_ = 0 ;
999  delete [] addedCuts_ ;
1000  addedCuts_ = NULL ;
1001/*
1002  Set up an empty heap and associated data structures to hold the live set
1003  (problems which require further exploration).
1004*/
1005  tree_->setComparison(*nodeCompare_) ;
1006/*
1007  Used to record the path from a node to the root of the search tree, so that
1008  we can then traverse from the root to the node when restoring a subproblem.
1009*/
1010  maximumDepth_ = 10 ;
1011  delete [] walkback_ ;
1012  walkback_ = new CbcNodeInfo * [maximumDepth_] ;
1013/*
1014  Used to generate bound edits for CbcPartialNodeInfo.
1015*/
1016  double * lowerBefore = new double [numberColumns] ;
1017  double * upperBefore = new double [numberColumns] ;
1018/*
1019 
1020  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
1021  lifting. It will iterate a generate/reoptimise loop (including reduced cost
1022  fixing) until no cuts are generated, the change in objective falls off,  or
1023  the limit on the number of rounds of cut generation is exceeded.
1024
1025  At the end of all this, any cuts will be recorded in cuts and also
1026  installed in the solver's constraint system. We'll have reoptimised, and
1027  removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
1028  adjusted accordingly).
1029
1030  Tell cut generators they can be a bit more aggressive at root node
1031
1032  TODO: Why don't we make a copy of the solution after solveWithCuts?
1033  TODO: If numberUnsatisfied == 0, don't we have a solution?
1034*/
1035  phase_=1;
1036  int iCutGenerator;
1037  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
1038    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
1039    generator->setAggressiveness(generator->getAggressiveness()+100);
1040  }
1041  OsiCuts cuts ;
1042  int anyAction = -1 ;
1043  numberOldActiveCuts_ = 0 ;
1044  numberNewCuts_ = 0 ;
1045  // Array to mark solution
1046  delete [] usedInSolution_;
1047  usedInSolution_ = new int[numberColumns];
1048  CoinZeroN(usedInSolution_,numberColumns);
1049/*
1050  For printing totals and for CbcNode (numberNodes_)
1051*/
1052  numberIterations_ = 0 ;
1053  numberNodes_ = 0 ;
1054  numberNodes2_ = 0 ;
1055  maximumStatistics_=0;
1056  maximumDepthActual_=0;
1057  numberDJFixed_=0.0;
1058  // Do heuristics
1059  {
1060    double * newSolution = new double [numberColumns] ;
1061    double heuristicValue = getCutoff() ;
1062    int found = -1; // no solution found
1063
1064    currentPassNumber_ = 1; // so root heuristics will run
1065    int i;
1066    for (i = 0;i<numberHeuristics_;i++) {
1067      // see if heuristic will do anything
1068      double saveValue = heuristicValue ;
1069      int ifSol = heuristic_[i]->solution(heuristicValue,
1070                                          newSolution);
1071      if (ifSol>0) {
1072        // better solution found
1073        found = i ;
1074        incrementUsed(newSolution);
1075        // increment number of solutions so other heuristics can test
1076        numberSolutions_++;
1077        numberHeuristicSolutions_++;
1078      } else {
1079        heuristicValue = saveValue ;
1080      }
1081    }
1082    currentPassNumber_ = 0;
1083    /*
1084      Did any of the heuristics turn up a new solution? Record it before we free
1085      the vector.
1086    */
1087    if (found >= 0) { 
1088      lastHeuristic_ = heuristic_[found];
1089      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
1090      CbcTreeLocal * tree
1091          = dynamic_cast<CbcTreeLocal *> (tree_);
1092      if (tree)
1093        tree->passInSolution(bestSolution_,heuristicValue);
1094    }
1095    for (i = 0;i<numberHeuristics_;i++) {
1096      // delete FPump
1097      CbcHeuristicFPump * pump
1098        = dynamic_cast<CbcHeuristicFPump *> (heuristic_[i]);
1099      if (pump) {
1100        delete pump;
1101        numberHeuristics_ --;
1102        for (int j=i;j<numberHeuristics_;j++)
1103          heuristic_[j] = heuristic_[j+1];
1104      }
1105    }
1106    delete [] newSolution ;
1107  } 
1108  statistics_ = NULL;
1109  // Do on switch
1110  if (doStatistics) {
1111    maximumStatistics_=10000;
1112    statistics_ = new CbcStatistics * [maximumStatistics_];
1113    memset(statistics_,0,maximumStatistics_*sizeof(CbcStatistics *));
1114  }
1115
1116  { int iObject ;
1117    int preferredWay ;
1118    int numberUnsatisfied = 0 ;
1119    memcpy(currentSolution_,solver_->getColSolution(),
1120           numberColumns*sizeof(double)) ;
1121    // point to useful information
1122    OsiBranchingInformation usefulInfo=usefulInformation();
1123
1124    for (iObject = 0 ; iObject < numberObjects_ ; iObject++)
1125    { double infeasibility =
1126          object_[iObject]->infeasibility(&usefulInfo,preferredWay) ;
1127      if (infeasibility ) numberUnsatisfied++ ; }
1128    // replace solverType
1129    if(solverCharacteristics_->tryCuts())  {
1130
1131      if (numberUnsatisfied)   {
1132        feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
1133                                 NULL);
1134      } else if (solverCharacteristics_->solutionAddsCuts()||
1135                 solverCharacteristics_->alwaysTryCutsAtRootNode()) {
1136        // may generate cuts and turn the solution
1137        //to an infeasible one
1138        feasible = solveWithCuts(cuts, 1,
1139                                 NULL);
1140      }
1141    }
1142    // check extra info on feasibility
1143    if (!solverCharacteristics_->mipFeasible())
1144      feasible = false;
1145  }
1146  // make cut generators less aggressive
1147  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
1148    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
1149    generator->setAggressiveness(generator->getAggressiveness()-100);
1150  }
1151  currentNumberCuts_ = numberNewCuts_ ;
1152  // See if can stop on gap
1153  stoppedOnGap_ = false ;
1154  bestPossibleObjective_ = solver_->getObjValue()*solver_->getObjSense();
1155  double testGap = CoinMax(dblParam_[CbcAllowableGap],
1156                           CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_))
1157                           *dblParam_[CbcAllowableFractionGap]);
1158  if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
1159    if (bestPossibleObjective_<getCutoff())
1160      stoppedOnGap_ = true ;
1161    feasible = false;
1162  }
1163/*
1164  We've taken the continuous relaxation as far as we can. Time to branch.
1165  The first order of business is to actually create a node. chooseBranch
1166  currently uses strong branching to evaluate branch object candidates,
1167  unless forced back to simple branching. If chooseBranch concludes that a
1168  branching candidate is monotone (anyAction == -1) or infeasible (anyAction
1169  == -2) when forced to integer values, it returns here immediately.
1170
1171  Monotone variables trigger a call to resolve(). If the problem remains
1172  feasible, try again to choose a branching variable. At the end of the loop,
1173  resolved == true indicates that some variables were fixed.
1174
1175  Loss of feasibility will result in the deletion of newNode.
1176*/
1177
1178  bool resolved = false ;
1179  CbcNode *newNode = NULL ;
1180  numberFixedAtRoot_=0;
1181  numberFixedNow_=0;
1182  int numberIterationsAtContinuous = numberIterations_;
1183  //solverCharacteristics_->setSolver(solver_);
1184  if (feasible) {
1185    newNode = new CbcNode ;
1186    // Set objective value (not so obvious if NLP etc)
1187    setObjectiveValue(newNode,NULL);
1188    anyAction = -1 ;
1189    // To make depth available we may need a fake node
1190    CbcNode fakeNode;
1191    if (!currentNode_) {
1192      // Not true if sub trees assert (!numberNodes_);
1193      currentNode_=&fakeNode;
1194    }
1195    phase_=3;
1196    // only allow 1000 passes
1197    int numberPassesLeft=1000;
1198    // This is first crude step
1199    if (numberAnalyzeIterations_) {
1200      delete [] analyzeResults_;
1201      analyzeResults_ = new double [4*numberIntegers_];
1202      numberFixedAtRoot_=newNode->analyze(this,analyzeResults_);
1203      if (numberFixedAtRoot_>0) {
1204        printf("%d fixed by analysis\n",numberFixedAtRoot_);
1205        setPointers(solver_);
1206        numberFixedNow_ = numberFixedAtRoot_;
1207      } else if (numberFixedAtRoot_<0) {
1208        printf("analysis found to be infeasible\n");
1209        anyAction=-2;
1210        delete newNode ;
1211        newNode = NULL ;
1212        feasible = false ;
1213      }
1214    }
1215    OsiSolverBranch * branches = NULL;
1216    anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved,
1217                             NULL,NULL,NULL,branches);
1218    if (anyAction == -2||newNode->objectiveValue() >= cutoff) {
1219      if (anyAction != -2) {
1220        // zap parent nodeInfo
1221#ifdef COIN_DEVELOP
1222        printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
1223#endif
1224        if (newNode->nodeInfo())
1225          newNode->nodeInfo()->nullParent();
1226      }
1227      delete newNode ;
1228      newNode = NULL ;
1229      feasible = false ;
1230    }
1231  }
1232/*
1233  At this point, the root subproblem is infeasible or fathomed by bound
1234  (newNode == NULL), or we're live with an objective value that satisfies the
1235  current objective cutoff.
1236*/
1237  assert (!newNode || newNode->objectiveValue() <= cutoff) ;
1238  // Save address of root node as we don't want to delete it
1239  // initialize for print out
1240  int lastDepth=0;
1241  int lastUnsatisfied=0;
1242  if (newNode)
1243    lastUnsatisfied=newNode->numberUnsatisfied();
1244/*
1245  The common case is that the lp relaxation is feasible but doesn't satisfy
1246  integrality (i.e., newNode->branchingObject(), indicating we've been able to
1247  select a branching variable). Remove any cuts that have gone slack due to
1248  forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
1249  basis (via createInfo()) and stash the new cuts in the nodeInfo (via
1250  addCuts()). If, by some miracle, we have an integral solution at the root
1251  (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
1252  a valid solution for use by setBestSolution().
1253*/
1254  CoinWarmStartBasis *lastws = NULL ;
1255  if (feasible && newNode->branchingObject())
1256  { if (resolved)
1257    { takeOffCuts(cuts,false,NULL) ;
1258#     ifdef CHECK_CUT_COUNTS
1259      { printf("Number of rows after chooseBranch fix (root)"
1260               "(active only) %d\n",
1261                numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ;
1262        const CoinWarmStartBasis* debugws =
1263          dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
1264        debugws->print() ;
1265        delete debugws ; }
1266#     endif
1267    }
1268  //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
1269    newNode->nodeInfo()->addCuts(cuts,
1270                                 newNode->numberBranches(),whichGenerator_) ;
1271    if (lastws) delete lastws ;
1272    lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
1273  }
1274/*
1275  Continuous data to be used later
1276*/
1277  continuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
1278  continuousInfeasibilities_ = 0 ;
1279  if (newNode)
1280  { continuousObjective_ = newNode->objectiveValue() ;
1281    delete [] continuousSolution_;
1282    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
1283                                             numberColumns);
1284    continuousInfeasibilities_ = newNode->numberUnsatisfied() ; }
1285/*
1286  Bound may have changed so reset in objects
1287*/
1288  { int i ;
1289    for (i = 0;i < numberObjects_;i++)
1290      object_[i]->resetBounds(solver_) ; }
1291/*
1292  Feasible? Then we should have either a live node prepped for future
1293  expansion (indicated by variable() >= 0), or (miracle of miracles) an
1294  integral solution at the root node.
1295
1296  initializeInfo sets the reference counts in the nodeInfo object.  Since
1297  this node is still live, push it onto the heap that holds the live set.
1298*/
1299  double bestValue = 0.0 ;
1300  if (newNode) {
1301    bestValue = newNode->objectiveValue();
1302    if (newNode->branchingObject()) {
1303      newNode->initializeInfo() ;
1304      tree_->push(newNode) ;
1305      if (statistics_) {
1306        if (numberNodes2_==maximumStatistics_) {
1307          maximumStatistics_ = 2*maximumStatistics_;
1308          CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
1309          memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
1310          memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
1311          delete [] statistics_;
1312          statistics_=temp;
1313        }
1314        assert (!statistics_[numberNodes2_]);
1315        statistics_[numberNodes2_]=new CbcStatistics(newNode);
1316      }
1317      numberNodes2_++;
1318#     ifdef CHECK_NODE
1319      printf("Node %x on tree\n",newNode) ;
1320#     endif
1321    } else {
1322      // continuous is integer
1323      double objectiveValue = newNode->objectiveValue();
1324      setBestSolution(CBC_SOLUTION,objectiveValue,
1325                      solver_->getColSolution()) ;
1326      delete newNode ;
1327      newNode = NULL ;
1328    }
1329  }
1330
1331  if (printFrequency_ <= 0) {
1332    printFrequency_ = 1000 ;
1333    if (getNumCols() > 2000)
1334      printFrequency_ = 100 ;
1335  }
1336  /*
1337    It is possible that strong branching fixes one variable and then the code goes round
1338    again and again.  This can take too long.  So we need to warn user - just once.
1339  */
1340  numberLongStrong_=0;
1341  double totalTime = 0.0;
1342#ifdef CBC_THREAD
1343  CbcNode * createdNode=NULL;
1344  CbcModel ** threadModel = NULL;
1345  pthread_t * threadId = NULL;
1346  int * threadCount = NULL;
1347  pthread_mutex_t mutex;
1348  pthread_cond_t condition_main;
1349  pthread_mutex_t condition_mutex;
1350  pthread_mutex_t * mutex2 = NULL;
1351  pthread_cond_t * condition2 = NULL;
1352  threadStruct * threadInfo = NULL;
1353  bool locked=false;
1354  int threadStats[6];
1355  memset(threadStats,0,sizeof(threadStats));
1356  double timeWaiting=0.0;
1357  // For now just one model
1358  if (numberThreads_) {
1359    threadId = new pthread_t [numberThreads_];
1360    threadCount = new int [numberThreads_];
1361    CoinZeroN(threadCount,numberThreads_);
1362    pthread_mutex_init(&mutex,NULL);
1363    pthread_cond_init(&condition_main,NULL);
1364    pthread_mutex_init(&condition_mutex,NULL);
1365    threadModel = new CbcModel * [numberThreads_+1];
1366    threadInfo = new threadStruct [numberThreads_+1];
1367    mutex2 = new pthread_mutex_t [numberThreads_];
1368    condition2 = new pthread_cond_t [numberThreads_];
1369    // we don't want a strategy object
1370    CbcStrategy * saveStrategy = strategy_;
1371    strategy_ = NULL;
1372    for (int i=0;i<numberThreads_;i++) {
1373      pthread_mutex_init(mutex2+i,NULL);
1374      pthread_cond_init(condition2+i,NULL);
1375      threadId[i]=0;
1376      threadInfo[i].baseModel=this;
1377      threadModel[i]=new CbcModel(*this);
1378#ifdef COIN_HAS_CLP
1379      // Solver may need to know about model
1380      CbcModel * thisModel = threadModel[i];
1381      CbcOsiSolver * solver =
1382              dynamic_cast<CbcOsiSolver *>(thisModel->solver()) ;
1383      if (solver)
1384        solver->setCbcModel(thisModel);
1385#endif
1386      mutex_ = (void *) (threadInfo+i);
1387      threadModel[i]->moveToModel(this,-1);
1388      threadInfo[i].thisModel=threadModel[i];
1389      threadInfo[i].node=NULL;
1390      threadInfo[i].createdNode=NULL;
1391      threadInfo[i].threadIdOfBase=pthread_self();
1392      threadInfo[i].mutex=&mutex;
1393      threadInfo[i].mutex2=mutex2+i;
1394      threadInfo[i].condition2=condition2+i;
1395      threadInfo[i].returnCode=-1;
1396      threadInfo[i].timeLocked=0.0;
1397      threadInfo[i].timeWaitingToLock=0.0;
1398      threadInfo[i].timeWaitingToStart=0.0;
1399      threadInfo[i].timeInThread=0.0;
1400      threadInfo[i].numberTimesLocked=0;
1401      threadInfo[i].numberTimesUnlocked=0;
1402      threadInfo[i].numberTimesWaitingToStart=0;
1403      threadInfo[i].locked=false;
1404#if CBC_THREAD_DEBUG
1405      threadInfo[i].threadNumber=i+2;
1406#endif
1407      pthread_create(threadId+i,NULL,doNodesThread,threadInfo+i);
1408    }
1409    strategy_ = saveStrategy;
1410    // Do a partial one for base model
1411    threadInfo[numberThreads_].baseModel=this;
1412    threadModel[numberThreads_]=this;
1413    mutex_ = (void *) (threadInfo+numberThreads_);
1414    threadInfo[numberThreads_].node=NULL;
1415    threadInfo[numberThreads_].mutex=&mutex;
1416    threadInfo[numberThreads_].condition2=&condition_main;
1417    threadInfo[numberThreads_].mutex2=&condition_mutex;
1418    threadInfo[numberThreads_].timeLocked=0.0;
1419    threadInfo[numberThreads_].timeWaitingToLock=0.0;
1420    threadInfo[numberThreads_].numberTimesLocked=0;
1421    threadInfo[numberThreads_].numberTimesUnlocked=0;
1422    threadInfo[numberThreads_].locked=false;
1423#if CBC_THREAD_DEBUG
1424    threadInfo[numberThreads_].threadNumber=1;
1425#endif
1426  }
1427#endif
1428/*
1429  At last, the actual branch-and-cut search loop, which will iterate until
1430  the live set is empty or we hit some limit (integrality gap, time, node
1431  count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
1432  solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
1433  add the node to the live set.
1434
1435  The first action is to winnow the live set to remove nodes which are worse
1436  than the current objective cutoff.
1437*/
1438  while (true) {
1439#ifdef CBC_THREAD
1440    if (!locked) {
1441      lockThread();
1442      locked=true;
1443    }
1444#endif
1445    if (tree_->empty()) {
1446#ifdef CBC_THREAD
1447      if (numberThreads_) {
1448        // may still be outstanding nodes
1449        int iThread;
1450        for (iThread=0;iThread<numberThreads_;iThread++) {
1451          if (threadId[iThread]) {
1452            if (threadInfo[iThread].returnCode==0) 
1453              break;
1454          }
1455        }
1456        if (iThread<numberThreads_) {
1457          unlockThread();
1458          locked = false;
1459          pthread_cond_signal(threadInfo[iThread].condition2); // unlock in case
1460#define CONDITION_WAIT
1461#ifdef CONDITION_WAIT
1462          while (true) {
1463            pthread_mutex_lock(&condition_mutex);
1464            struct timespec absTime;
1465            clock_gettime(CLOCK_REALTIME,&absTime);
1466            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
1467            absTime.tv_nsec += 1000000; // millisecond
1468            if (absTime.tv_nsec>=1000000000) {
1469              absTime.tv_nsec -= 1000000000;
1470              absTime.tv_sec++;
1471            }
1472            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
1473            clock_gettime(CLOCK_REALTIME,&absTime);
1474            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
1475            timeWaiting += time2-time;
1476            pthread_mutex_unlock(&condition_mutex);
1477            if (threadInfo[iThread].returnCode!=0) 
1478              break;
1479          }
1480#else
1481          while (threadInfo[iThread].returnCode==0) {
1482            struct timespec moreTime;
1483            moreTime.tv_nsec = 10000;
1484            moreTime.tv_sec = 0;
1485            nanosleep(&moreTime,NULL);
1486          }
1487#endif
1488          threadModel[iThread]->moveToModel(this,1);
1489          assert (threadInfo[iThread].returnCode==1);
1490          // say available
1491          threadInfo[iThread].returnCode=-1;
1492          threadStats[4]++;
1493          continue;
1494        }
1495      }
1496      unlockThread();
1497      locked=false; // not needed as break
1498#endif
1499      break;
1500    }
1501#ifdef CBC_THREAD
1502    unlockThread();
1503    locked = false;
1504#endif
1505/*
1506  Check for abort on limits: node count, solution count, time, integrality gap.
1507*/
1508    totalTime = CoinCpuTime()-dblParam_[CbcStartSeconds] ;
1509    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
1510          numberSolutions_ < intParam_[CbcMaxNumSol] &&
1511          totalTime < dblParam_[CbcMaximumSeconds] &&
1512          !stoppedOnGap_&&!eventHappened_)) {
1513      // out of loop
1514      break;
1515    }
1516#ifdef BONMIN
1517    assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
1518#endif
1519    if (cutoff > getCutoff()) {
1520      double newCutoff = getCutoff();
1521      if (analyzeResults_) {
1522        // see if we could fix any (more)
1523        int n=0;
1524        double * newLower = analyzeResults_;
1525        double * objLower = newLower+numberIntegers_;
1526        double * newUpper = objLower+numberIntegers_;
1527        double * objUpper = newUpper+numberIntegers_;
1528        for (int i=0;i<numberIntegers_;i++) {
1529          if (objLower[i]>newCutoff) {
1530            n++;
1531            if (objUpper[i]>newCutoff) {
1532              newCutoff = -COIN_DBL_MAX;
1533              break;
1534            }
1535          } else if (objUpper[i]>newCutoff) {
1536            n++;
1537          }
1538        }
1539        if (newCutoff==-COIN_DBL_MAX) {
1540          printf("Root analysis says finished\n");
1541        } else if (n>numberFixedNow_) {
1542          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
1543          numberFixedNow_=n;
1544        }
1545      }
1546      if (eventHandler) {
1547        if (!eventHandler->event(CbcEventHandler::solution)) {
1548          eventHappened_=true; // exit
1549        }
1550      }
1551      lockThread();
1552      // Do from deepest
1553      tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ;
1554      nodeCompare_->newSolution(this) ;
1555      nodeCompare_->newSolution(this,continuousObjective_,
1556                                continuousInfeasibilities_) ;
1557      tree_->setComparison(*nodeCompare_) ;
1558      if (tree_->empty()) {
1559        unlockThread();
1560        break; // finished
1561      }
1562      unlockThread();
1563    }
1564    cutoff = getCutoff() ;
1565/*
1566    Periodic activities: Opportunities to
1567    + tweak the nodeCompare criteria,
1568    + check if we've closed the integrality gap enough to quit,
1569    + print a summary line to let the user know we're working
1570*/
1571    if ((numberNodes_%1000) == 0) {
1572      lockThread();
1573      bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ;
1574#ifdef CHECK_CUT_SIZE
1575      verifyCutSize (tree_, *this);
1576#endif
1577      // redo tree if wanted
1578      if (redoTree)
1579        tree_->setComparison(*nodeCompare_) ;
1580      unlockThread();
1581    }
1582    if (saveCompare&&!hotstartSolution_) {
1583      // hotstart switched off
1584      delete nodeCompare_; // off depth first
1585      nodeCompare_=saveCompare;
1586      saveCompare=NULL;
1587      // redo tree
1588      lockThread();
1589      tree_->setComparison(*nodeCompare_) ;
1590      unlockThread();
1591    }
1592    if ((numberNodes_%printFrequency_) == 0) {
1593      lockThread();
1594      int j ;
1595      int nNodes = tree_->size() ;
1596      bestPossibleObjective_ = 1.0e100 ;
1597      for (j = 0;j < nNodes;j++) {
1598        CbcNode * node = tree_->nodePointer(j) ;
1599        if (node&&node->objectiveValue() < bestPossibleObjective_)
1600          bestPossibleObjective_ = node->objectiveValue() ;
1601      }
1602      unlockThread();
1603      if (!intParam_[CbcPrinting]) {
1604        messageHandler()->message(CBC_STATUS,messages())
1605          << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
1606          <<getCurrentSeconds()
1607          << CoinMessageEol ;
1608      } else {
1609        messageHandler()->message(CBC_STATUS2,messages())
1610          << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
1611          <<lastDepth<<lastUnsatisfied<<numberIterations_
1612          <<getCurrentSeconds()
1613          << CoinMessageEol ;
1614      }
1615      if (!eventHandler->event(CbcEventHandler::treeStatus)) {
1616        eventHappened_=true; // exit
1617      }
1618    }
1619    // If no solution but many nodes - signal change in strategy
1620    if (numberNodes_>2*numberObjects_+1000&&stateOfSearch_!=2)
1621      stateOfSearch_=3;
1622    // See if can stop on gap
1623    double testGap = CoinMax(dblParam_[CbcAllowableGap],
1624                             CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_))
1625                             *dblParam_[CbcAllowableFractionGap]);
1626    if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
1627      stoppedOnGap_ = true ;
1628    }
1629   
1630#   ifdef CHECK_NODE_FULL
1631    verifyTreeNodes(tree_,*this) ;
1632#   endif
1633#   ifdef CHECK_CUT_COUNTS
1634    verifyCutCounts(tree_,*this) ;
1635#   endif
1636/*
1637  Now we come to the meat of the loop. To create the active subproblem, we'll
1638  pop the most promising node in the live set, rebuild the subproblem it
1639  represents, and then execute the current arm of the branch to create the
1640  active subproblem.
1641*/
1642    CbcNode *node = tree_->bestNode(cutoff) ;
1643    // Possible one on tree worse than cutoff
1644    if (!node)
1645      continue;
1646#ifndef CBC_THREAD
1647    int currentNumberCuts = 0 ;
1648    currentNode_=node; // so can be accessed elsewhere
1649#ifdef CBC_DEBUG
1650    printf("%d unsat, way %d, obj %g est %g\n",
1651           node->numberUnsatisfied(),node->way(),node->objectiveValue(),
1652           node->guessedObjectiveValue());
1653#endif
1654    // Save clone in branching decision
1655    if(branchingMethod_)
1656      branchingMethod_->saveBranchingObject(node->modifiableBranchingObject());
1657    // Say not on optimal path
1658    bool onOptimalPath=false;
1659#   ifdef CHECK_NODE
1660    printf("Node %x popped from tree - %d left, %d count\n",node,
1661           node->nodeInfo()->numberBranchesLeft(),
1662           node->nodeInfo()->numberPointingToThis()) ;
1663    printf("\tdepth = %d, z =  %g, unsat = %d, var = %d.\n",
1664           node->depth(),node->objectiveValue(),
1665           node->numberUnsatisfied(),
1666           node->columnNumber()) ;
1667#   endif
1668    lastDepth=node->depth();
1669    lastUnsatisfied=node->numberUnsatisfied();
1670
1671/*
1672  Rebuild the subproblem for this node:  Call addCuts() to adjust the model
1673  to recreate the subproblem for this node (set proper variable bounds, add
1674  cuts, create a basis).  This may result in the problem being fathomed by
1675  bound or infeasibility. Returns 1 if node is fathomed.
1676  Execute the current arm of the branch: If the problem survives, save the
1677  resulting variable bounds and call branch() to modify variable bounds
1678  according to the current arm of the branching object. If we're processing
1679  the final arm of the branching object, flag the node for removal from the
1680  live set.
1681*/
1682    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
1683    newNode = NULL ;
1684    int branchesLeft=0;
1685    if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_))
1686    { int i ;
1687      const double * lower = getColLower() ;
1688      const double * upper = getColUpper() ;
1689      for (i = 0 ; i < numberColumns ; i++)
1690      { lowerBefore[i]= lower[i] ;
1691        upperBefore[i]= upper[i] ; }
1692      if ((solverCharacteristics_->extraCharacteristics()&2)!=0) {
1693        solverCharacteristics_->setBeforeLower(lowerBefore);
1694        solverCharacteristics_->setBeforeUpper(upperBefore);
1695      }
1696      if (messageHandler()->logLevel()>2)
1697        node->modifiableBranchingObject()->print();
1698      if (!useOsiBranching) 
1699        branchesLeft = node->branch(NULL); // old way
1700      else
1701        branchesLeft = node->branch(solver_); // new way
1702      if (branchesLeft) {
1703        // set nodenumber correctly
1704        node->nodeInfo()->setNodeNumber(numberNodes2_);
1705        tree_->push(node) ;
1706        if (statistics_) {
1707          if (numberNodes2_==maximumStatistics_) {
1708            maximumStatistics_ = 2*maximumStatistics_;
1709            CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
1710            memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
1711            memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
1712            delete [] statistics_;
1713            statistics_=temp;
1714          }
1715          assert (!statistics_[numberNodes2_]);
1716          statistics_[numberNodes2_]=new CbcStatistics(node);
1717        }
1718        numberNodes2_++;
1719        //nodeOnTree=true; // back on tree
1720        //deleteNode = false ;
1721#       ifdef CHECK_NODE
1722        printf("Node %x pushed back on tree - %d left, %d count\n",node,
1723               nodeInfo->numberBranchesLeft(),
1724               nodeInfo->numberPointingToThis()) ;
1725#       endif
1726      } else {
1727        //deleteNode = true ;
1728        if (!nodeInfo->numberBranchesLeft())
1729          nodeInfo->allBranchesGone(); // can clean up
1730      }
1731      if ((specialOptions_&1)!=0) {
1732        /*
1733          This doesn't work as intended --- getRowCutDebugger will return null
1734          unless the current feasible solution region includes the optimal solution
1735          that RowCutDebugger knows. There's no way to tell inactive from off the
1736          optimal path.
1737        */
1738        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
1739        if (debugger) {
1740          onOptimalPath=true;
1741          printf("On optimal path\n") ;
1742        }
1743      }
1744     
1745/*
1746  Reoptimize, possibly generating cuts and/or using heuristics to find
1747  solutions.  Cut reference counts are unaffected unless we lose feasibility,
1748  in which case solveWithCuts() will make the adjustment.
1749*/
1750      phase_=2;
1751      cuts = OsiCuts() ;
1752      currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_ ;
1753      int saveNumber = numberIterations_;
1754      if(solverCharacteristics_->solutionAddsCuts()) {
1755        int returnCode=resolve(node ? node->nodeInfo() : NULL,1);
1756        feasible = returnCode != 0;
1757        if (feasible) {
1758          int iObject ;
1759          int preferredWay ;
1760          int numberUnsatisfied = 0 ;
1761          memcpy(currentSolution_,solver_->getColSolution(),
1762                 numberColumns*sizeof(double)) ;
1763          // point to useful information
1764          OsiBranchingInformation usefulInfo=usefulInformation();
1765         
1766          for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1767            double infeasibility =
1768              object_[iObject]->infeasibility(&usefulInfo,preferredWay) ;
1769            if (infeasibility ) numberUnsatisfied++ ;
1770          }
1771          if (returnCode>0) {
1772            if (numberUnsatisfied)   {
1773              feasible = solveWithCuts(cuts,maximumCutPasses_,node);
1774            } else {
1775              // may generate cuts and turn the solution
1776              //to an infeasible one
1777              feasible = solveWithCuts(cuts, 1,
1778                                       node);
1779#if 0
1780              currentNumberCuts_ = cuts.sizeRowCuts();
1781              if (currentNumberCuts_ >= maximumNumberCuts_) {
1782                maximumNumberCuts_ = currentNumberCuts;
1783                delete [] addedCuts_;
1784                addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
1785              }
1786#endif
1787            }
1788          }
1789          // check extra info on feasibility
1790          if (!solverCharacteristics_->mipFeasible()) {
1791            feasible = false;
1792            solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
1793          }
1794        }
1795      } else {
1796        // normal
1797        //int zzzzzz=0;
1798        //if (zzzzzz)
1799        //solver_->writeMps("before");
1800        feasible = solveWithCuts(cuts,maximumCutPasses_,node);
1801      }
1802      if ((specialOptions_&1)!=0&&onOptimalPath) {
1803        if (!solver_->getRowCutDebugger()) {
1804          // dj fix did something???
1805          solver_->writeMps("infeas2");
1806          assert (solver_->getRowCutDebugger()) ;
1807        }
1808      }
1809      if (statistics_) {
1810        assert (numberNodes2_);
1811        assert (statistics_[numberNodes2_-1]);
1812        assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
1813        statistics_[numberNodes2_-1]->endOfBranch(numberIterations_-saveNumber,
1814                                               feasible ? solver_->getObjValue()
1815                                               : COIN_DBL_MAX);
1816      }
1817/*
1818  Are we still feasible? If so, create a node and do the work to attach a
1819  branching object, reoptimising as needed if chooseBranch() identifies
1820  monotone objects.
1821
1822  Finally, attach a partial nodeInfo object and store away any cuts that we
1823  created back in solveWithCuts. addCuts() will initialise the reference
1824  counts for these new cuts.
1825
1826  This next test can be problematic if we've discovered an
1827  alternate equivalent answer and subsequently fathom the solution
1828  known to the row cut debugger due to bounds.
1829*/
1830        if (onOptimalPath) {
1831          bool objLim = solver_->isDualObjectiveLimitReached() ;
1832          if (!feasible && !objLim) {
1833            printf("infeas2\n");
1834            solver_->writeMps("infeas");
1835            CoinWarmStartBasis *slack =
1836              dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
1837            solver_->setWarmStart(slack);
1838            delete slack ;
1839            solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
1840            solver_->initialSolve();
1841            assert (!solver_->isProvenOptimal());
1842          }
1843          assert (feasible || objLim);
1844        }
1845        bool checkingNode=false;
1846        if (feasible) {
1847          newNode = new CbcNode ;
1848          // Set objective value (not so obvious if NLP etc)
1849          setObjectiveValue(newNode,node);
1850          anyAction =-1 ;
1851          resolved = false ;
1852          if (newNode->objectiveValue() >= getCutoff()) 
1853            anyAction=-2;
1854          // only allow at most a few passes
1855          int numberPassesLeft=5;
1856          checkingNode=true;
1857        OsiSolverBranch * branches=NULL;
1858        // point to useful information
1859        anyAction = chooseBranch(newNode, numberPassesLeft,node, cuts,resolved,
1860                                 lastws, lowerBefore, upperBefore, branches);
1861/*
1862  If we end up infeasible, we can delete the new node immediately. Since this
1863  node won't be needing the cuts we collected, decrement the reference counts.
1864  If we are feasible, then we'll be placing this node into the live set, so
1865  increment the reference count in the current (parent) nodeInfo.
1866*/
1867        if (anyAction == -2)
1868          { delete newNode ;
1869          newNode = NULL ;
1870          // say strong doing well
1871          if (checkingNode)
1872            setSpecialOptions(specialOptions_|8);
1873          for (i = 0 ; i < currentNumberCuts_ ; i++)
1874            { if (addedCuts_[i])
1875              { if (!addedCuts_[i]->decrement(1))
1876                delete addedCuts_[i] ; } } }
1877        else
1878          { nodeInfo->increment() ;
1879          if ((numberNodes_%20)==0) {
1880            // say strong not doing as well
1881            setSpecialOptions(specialOptions_&~8);
1882          }
1883        }
1884        }
1885/*
1886  At this point, there are three possibilities:
1887    * newNode is live and will require further branching to resolve
1888      (variable() >= 0). Increment the cut reference counts by
1889      numberBranches() to allow for use by children of this node, and
1890      decrement by 1 because we've executed one arm of the branch of our
1891      parent (consuming one reference). Before we push newNode onto the
1892      search tree, try for a heuristic solution.
1893    * We have a solution, in which case newNode is non-null but we have no
1894      branching variable. Decrement the cut counts and save the solution.
1895    * The node was found to be infeasible, in which case it's already been
1896      deleted, and newNode is null.
1897*/
1898        if (!eventHandler->event(CbcEventHandler::node)) {
1899          eventHappened_=true; // exit
1900        }
1901        assert (!newNode || newNode->objectiveValue() <= getCutoff()) ;
1902        if (statistics_) {
1903          assert (numberNodes2_);
1904          assert (statistics_[numberNodes2_-1]);
1905          assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
1906          if (newNode)
1907            statistics_[numberNodes2_-1]->updateInfeasibility(newNode->numberUnsatisfied());
1908          else
1909            statistics_[numberNodes2_-1]->sayInfeasible();
1910        }
1911        if (newNode) {
1912          if (newNode->branchingObject() == NULL&&solverCharacteristics_->solverType()==4) {
1913            // need to check if any cuts would do anything
1914            OsiCuts theseCuts;
1915            // reset probing info
1916            if (probingInfo_)
1917              probingInfo_->initializeFixing();
1918            for (int i = 0;i<numberCutGenerators_;i++) {
1919              bool generate = generator_[i]->normal();
1920              // skip if not optimal and should be (maybe a cut generator has fixed variables)
1921              if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
1922                generate=false;
1923              if (!generator_[i]->mustCallAgain())
1924                generate=false; // only special cuts
1925              if (generate) {
1926                generator_[i]->generateCuts(theseCuts,true,solver_,NULL) ;
1927                int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
1928                if (numberRowCutsAfter) {
1929                  // need dummy branch
1930                  newNode->setBranchingObject(new CbcDummyBranchingObject(this));
1931                  newNode->nodeInfo()->initializeInfo(1);
1932                  break;
1933                }
1934              }
1935            }
1936          }
1937          if (newNode->branchingObject())
1938          { handler_->message(CBC_BRANCH,messages_)
1939               << numberNodes_<< newNode->objectiveValue()
1940               << newNode->numberUnsatisfied()<< newNode->depth()
1941               << CoinMessageEol ;
1942            // Increment cut counts (taking off current)
1943            int numberLeft = newNode->numberBranches() ;
1944            for (i = 0;i < currentNumberCuts_;i++)
1945            { if (addedCuts_[i])
1946              {
1947#               ifdef CHECK_CUT_COUNTS
1948                printf("Count on cut %x increased by %d\n",addedCuts_[i],
1949                        numberLeft-1) ;
1950#               endif
1951                addedCuts_[i]->increment(numberLeft-1) ; } }
1952
1953            double estValue = newNode->guessedObjectiveValue() ;
1954            int found = -1 ;
1955            // no - overhead on small problems solver_->resolve() ;     // double check current optimal
1956            // assert (!solver_->getIterationCount());
1957            double * newSolution = new double [numberColumns] ;
1958            double heurValue = getCutoff() ;
1959            int iHeur ;
1960            for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++)
1961            { double saveValue = heurValue ;
1962              int ifSol = heuristic_[iHeur]->solution(heurValue,newSolution) ;
1963              if (ifSol > 0) {
1964                // new solution found
1965                found = iHeur ;
1966                incrementUsed(newSolution);
1967              }
1968              else
1969              if (ifSol < 0)    // just returning an estimate
1970              { estValue = CoinMin(heurValue,estValue) ;
1971                heurValue = saveValue ; } }
1972            if (found >= 0) {
1973              lastHeuristic_ = heuristic_[found];
1974              setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
1975            }
1976            delete [] newSolution ;
1977            newNode->setGuessedObjectiveValue(estValue) ;
1978            tree_->push(newNode) ;
1979            if (statistics_) {
1980              if (numberNodes2_==maximumStatistics_) {
1981                maximumStatistics_ = 2*maximumStatistics_;
1982                CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
1983                memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
1984                memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
1985                delete [] statistics_;
1986                statistics_=temp;
1987              }
1988              assert (!statistics_[numberNodes2_]);
1989              statistics_[numberNodes2_]=new CbcStatistics(newNode);
1990            }
1991            numberNodes2_++;
1992#           ifdef CHECK_NODE
1993            printf("Node %x pushed on tree c\n",newNode) ;
1994#           endif
1995          }
1996          else
1997          { 
1998            if(solverCharacteristics_ && //we may be in a non standard bab
1999               solverCharacteristics_->solutionAddsCuts()// we are in some kind of OA based bab.
2000               )
2001              {
2002                std::cerr<<"You should never get here"<<std::endl;
2003                throw CoinError("Nodes should not be fathomed on integer infeasibility in this setting",
2004                                "branchAndBound","CbcModel") ;
2005              }
2006            for (i = 0 ; i < currentNumberCuts_ ; i++)
2007            { if (addedCuts_[i])
2008              { if (!addedCuts_[i]->decrement(1))
2009                  delete addedCuts_[i] ; } }
2010          double objectiveValue = newNode->objectiveValue();
2011            setBestSolution(CBC_SOLUTION,objectiveValue,
2012                            solver_->getColSolution()) ;
2013            lastHeuristic_ = NULL;
2014            incrementUsed(solver_->getColSolution());
2015            //assert(nodeInfo->numberPointingToThis() <= 2) ;
2016            // avoid accidental pruning, if newNode was final branch arm
2017            nodeInfo->increment();
2018            delete newNode ;
2019            nodeInfo->decrement() ; } }
2020/*
2021  This node has been completely expanded and can be removed from the live
2022  set.
2023*/
2024      if (branchesLeft)
2025      { 
2026      }
2027      else
2028      { 
2029        if (!nodeInfo->numberBranchesLeft())
2030          nodeInfo->allBranchesGone(); // can clean up
2031        delete node ; }
2032    } else {
2033      // add cuts found to be infeasible (on bound)!
2034      abort();
2035      delete node;
2036    }
2037/*
2038  Delete cuts to get back to the original system.
2039
2040  I'm thinking this is redundant --- the call to addCuts that conditions entry
2041  to this code block also performs this action.
2042*/
2043      int numberToDelete = getNumRows()-numberRowsAtContinuous_ ;
2044      if (numberToDelete)
2045      { int * delRows = new int[numberToDelete] ;
2046        int i ;
2047        for (i = 0 ; i < numberToDelete ; i++)
2048        { delRows[i] = i+numberRowsAtContinuous_ ; }
2049        solver_->deleteRows(numberToDelete,delRows) ;
2050        delete [] delRows ; }
2051#else // end of not CBC_THREAD
2052      if (!numberThreads_) {
2053        doOneNode(this,node,createdNode);
2054      } else {
2055        threadStats[0]++;
2056        //need to think
2057        int iThread;
2058        // Start one off if any available
2059        for (iThread=0;iThread<numberThreads_;iThread++) {
2060          if (threadInfo[iThread].returnCode==-1) 
2061            break;
2062        }
2063        if (iThread<numberThreads_) {
2064          threadInfo[iThread].node=node;
2065          assert (threadInfo[iThread].returnCode==-1);
2066          // say in use
2067          threadInfo[iThread].returnCode=0;
2068          threadModel[iThread]->moveToModel(this,0);
2069          pthread_cond_signal(threadInfo[iThread].condition2); // unlock
2070          threadCount[iThread]++;
2071        }
2072        lockThread();
2073        locked=true;
2074        // see if any finished
2075        for (iThread=0;iThread<numberThreads_;iThread++) {
2076          if (threadInfo[iThread].returnCode>0) 
2077            break;
2078        }
2079        unlockThread();
2080        locked=false;
2081        if (iThread<numberThreads_) {
2082          threadModel[iThread]->moveToModel(this,1);
2083          assert (threadInfo[iThread].returnCode==1);
2084          // say available
2085          threadInfo[iThread].returnCode=-1;
2086          // carry on
2087          threadStats[3]++;
2088        } else {
2089          // Start one off if any available
2090          for (iThread=0;iThread<numberThreads_;iThread++) {
2091            if (threadInfo[iThread].returnCode==-1) 
2092              break;
2093          }
2094          if (iThread<numberThreads_) {
2095            lockThread();
2096            locked=true;
2097            // If any on tree get
2098            if (!tree_->empty()) {
2099              //node = tree_->bestNode(cutoff) ;
2100              //assert (node);
2101              threadStats[1]++;
2102              continue; // ** get another node
2103            }
2104            unlockThread();
2105            locked=false;
2106          }
2107          // wait (for debug could sleep and use test)
2108          bool finished=false;
2109          while (!finished) {
2110#ifdef CONDITION_WAIT
2111            pthread_mutex_lock(&condition_mutex);
2112            struct timespec absTime;
2113            clock_gettime(CLOCK_REALTIME,&absTime);
2114            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
2115            absTime.tv_nsec += 1000000; // millisecond
2116            if (absTime.tv_nsec>=1000000000) {
2117              absTime.tv_nsec -= 1000000000;
2118              absTime.tv_sec++;
2119            }
2120            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
2121            clock_gettime(CLOCK_REALTIME,&absTime);
2122            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
2123            timeWaiting += time2-time;
2124            pthread_mutex_unlock(&condition_mutex);
2125#else
2126            struct timespec moreTime;
2127            moreTime.tv_nsec = 100;
2128            moreTime.tv_sec = 0;
2129            nanosleep(&moreTime,NULL);
2130#endif
2131            for (iThread=0;iThread<numberThreads_;iThread++) {
2132              if (threadInfo[iThread].returnCode>0) {
2133                finished=true;
2134                break;
2135              } else if (threadInfo[iThread].returnCode==0) {
2136                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
2137              }
2138            }
2139          }
2140          assert (iThread<numberThreads_);
2141          threadModel[iThread]->moveToModel(this,1);
2142          node = threadInfo[iThread].node;
2143          threadInfo[iThread].node=NULL;
2144          assert (threadInfo[iThread].returnCode==1);
2145          // say available
2146          threadInfo[iThread].returnCode=-1;
2147          // carry on
2148          threadStats[2]++;
2149        }
2150      }
2151      //lastDepth=node->depth();
2152      //lastUnsatisfied=node->numberUnsatisfied();
2153#endif // end of CBC_THREAD
2154  }
2155#ifdef CBC_THREAD
2156  if (numberThreads_) {
2157    //printf("stats ");
2158    //for (unsigned int j=0;j<sizeof(threadStats)/sizeof(int);j++)
2159    //printf("%d ",threadStats[j]);
2160    //printf("\n");
2161    int i;
2162    // Seems to be bug in CoinCpu on Linux - does threads as well despite documentation
2163    double time=0.0;
2164    for (i=0;i<numberThreads_;i++) 
2165      time += threadInfo[i].timeInThread;
2166    bool goodTimer = time<(CoinCpuTime() - dblParam_[CbcStartSeconds]);
2167    for (i=0;i<numberThreads_;i++) {
2168      pthread_cond_signal(threadInfo[i].condition2); // unlock
2169      while (threadInfo[i].returnCode==0) {
2170        struct timespec moreTime;
2171        moreTime.tv_nsec = 10000;
2172        moreTime.tv_sec = 0;
2173        nanosleep(&moreTime,NULL);
2174      }
2175      threadModel[i]->numberThreads_=0; // say exit
2176      threadInfo[i].returnCode=0;
2177      pthread_cond_signal(threadInfo[i].condition2); // unlock
2178      pthread_join(threadId[i],NULL);
2179      threadModel[i]->moveToModel(this,2);
2180      pthread_mutex_destroy (threadInfo[i].mutex2);
2181      pthread_cond_destroy (threadInfo[i].condition2);
2182      assert (threadInfo[i].numberTimesLocked==threadInfo[i].numberTimesUnlocked);
2183      handler_->message(CBC_THREAD_STATS,messages_)
2184        <<"Thread";
2185      handler_->printing(true)
2186        <<i<<threadCount[i]<<threadInfo[i].timeWaitingToStart;
2187      handler_->printing(goodTimer)<<threadInfo[i].timeInThread;
2188      handler_->printing(false)<<0.0;
2189      handler_->printing(true)<<threadInfo[i].numberTimesLocked
2190        <<threadInfo[i].timeLocked<<threadInfo[i].timeWaitingToLock
2191        <<CoinMessageEol;
2192    }
2193    assert (threadInfo[numberThreads_].numberTimesLocked==threadInfo[numberThreads_].numberTimesUnlocked);
2194    handler_->message(CBC_THREAD_STATS,messages_)
2195      <<"Main thread";
2196    handler_->printing(false)<<0<<0<<0.0;
2197    handler_->printing(false)<<0.0;
2198    handler_->printing(true)<<timeWaiting;
2199    handler_->printing(true)<<threadInfo[numberThreads_].numberTimesLocked
2200      <<threadInfo[numberThreads_].timeLocked<<threadInfo[numberThreads_].timeWaitingToLock
2201      <<CoinMessageEol;
2202    pthread_mutex_destroy (&mutex);
2203    pthread_cond_destroy (&condition_main);
2204    pthread_mutex_destroy (&condition_mutex);
2205    // delete models (here in case some point to others)
2206    for (i=0;i<numberThreads_;i++) {
2207      delete threadModel[i];
2208    }
2209    delete [] mutex2;
2210    delete [] condition2;
2211    delete [] threadId;
2212    delete [] threadInfo;
2213    delete [] threadModel;
2214    delete [] threadCount;
2215    mutex_=NULL;
2216  }
2217#endif
2218/*
2219  End of the non-abort actions. The next block of code is executed if we've
2220  aborted because we hit one of the limits. Clean up by deleting the live set
2221  and break out of the node processing loop. Note that on an abort, node may
2222  have been pushed back onto the tree for further processing, in which case
2223  it'll be deleted in cleanTree. We need to check.
2224*/
2225    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
2226        numberSolutions_ < intParam_[CbcMaxNumSol] &&
2227        totalTime < dblParam_[CbcMaximumSeconds] &&
2228        !stoppedOnGap_&&!eventHappened_)) {
2229      if (tree_->size())
2230        tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
2231      delete nextRowCut_;
2232      if (stoppedOnGap_)
2233        { messageHandler()->message(CBC_GAP,messages())
2234          << bestObjective_-bestPossibleObjective_
2235          << dblParam_[CbcAllowableGap]
2236          << dblParam_[CbcAllowableFractionGap]*100.0
2237          << CoinMessageEol ;
2238        secondaryStatus_ = 2;
2239        status_ = 0 ; }
2240        else
2241          if (isNodeLimitReached())
2242            { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
2243            secondaryStatus_ = 3;
2244            status_ = 1 ; }
2245          else
2246        if (totalTime >= dblParam_[CbcMaximumSeconds])
2247          { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
2248          secondaryStatus_ = 4;
2249          status_ = 1 ; }
2250        else
2251          if (eventHappened_)
2252            { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
2253            secondaryStatus_ = 5;
2254            status_ = 5 ; }
2255          else
2256            { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
2257            secondaryStatus_ = 6;
2258            status_ = 1 ; }
2259    }
2260/*
2261  That's it, we've exhausted the search tree, or broken out of the loop because
2262  we hit some limit on evaluation.
2263
2264  We may have got an intelligent tree so give it one more chance
2265*/
2266  // Tell solver we are not in Branch and Cut
2267  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
2268  tree_->endSearch();
2269  //  If we did any sub trees - did we give up on any?
2270  if ( numberStoppedSubTrees_)
2271    status_=1;
2272  if (!status_) {
2273    // Set best possible unless stopped on gap
2274    if(secondaryStatus_ != 2)
2275      bestPossibleObjective_=bestObjective_;
2276    handler_->message(CBC_END_GOOD,messages_)
2277      << bestObjective_ << numberIterations_ << numberNodes_<<getCurrentSeconds()
2278      << CoinMessageEol ;
2279  } else {
2280    handler_->message(CBC_END,messages_)
2281      << bestObjective_ <<bestPossibleObjective_
2282      << numberIterations_ << numberNodes_<<getCurrentSeconds()
2283      << CoinMessageEol ;
2284  }
2285  if (numberStrongIterations_)
2286    handler_->message(CBC_STRONG_STATS,messages_)
2287      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
2288      << strongInfo_[1] << CoinMessageEol ;
2289  handler_->message(CBC_OTHER_STATS,messages_)
2290      << maximumDepthActual_
2291      << numberDJFixed_ << CoinMessageEol ;
2292  if (statistics_) {
2293    // report in some way
2294    int * lookup = new int[numberObjects_];
2295    int i;
2296    for (i=0;i<numberObjects_;i++) 
2297      lookup[i]=-1;
2298    bool goodIds=true;
2299    for (i=0;i<numberObjects_;i++) {
2300      int iColumn = object_[i]->columnNumber();
2301      if(iColumn>=0&&iColumn<numberColumns) {
2302        if (lookup[i]==-1) {
2303          lookup[i]=iColumn;
2304        } else {
2305          goodIds=false;
2306          break;
2307        }
2308      } else {
2309        goodIds=false;
2310        break;
2311      }
2312    }
2313    if (!goodIds) {
2314      delete [] lookup;
2315      lookup=NULL;
2316    }
2317    if (doStatistics==3) {
2318      printf("  node parent depth column   value                    obj      inf\n");
2319      for ( i=0;i<numberNodes2_;i++) {
2320        statistics_[i]->print(lookup);
2321      }
2322    }
2323    if (doStatistics>1) {
2324      // Find last solution
2325      int k;
2326      for (k=numberNodes2_-1;k>=0;k--) {
2327        if (statistics_[k]->endingObjective()!=COIN_DBL_MAX&&
2328            !statistics_[k]->endingInfeasibility())
2329          break;
2330      }
2331      if (k>=0) {
2332        int depth=statistics_[k]->depth();
2333        int * which = new int[depth+1];
2334        for (i=depth;i>=0;i--) {
2335          which[i]=k;
2336          k=statistics_[k]->parentNode();
2337        }
2338        printf("  node parent depth column   value                    obj      inf\n");
2339        for (i=0;i<=depth;i++) {
2340          statistics_[which[i]]->print(lookup);
2341        }
2342        delete [] which;
2343      }
2344    }
2345    // now summary
2346    int maxDepth=0;
2347    double averageSolutionDepth=0.0;
2348    int numberSolutions=0;
2349    double averageCutoffDepth=0.0;
2350    double averageSolvedDepth=0.0;
2351    int numberCutoff=0;
2352    int numberDown=0;
2353    int numberFirstDown=0;
2354    double averageInfDown=0.0;
2355    double averageObjDown=0.0;
2356    int numberCutoffDown=0;
2357    int numberUp=0;
2358    int numberFirstUp=0;
2359    double averageInfUp=0.0;
2360    double averageObjUp=0.0;
2361    int numberCutoffUp=0;
2362    double averageNumberIterations1=0.0;
2363    double averageValue=0.0;
2364    for ( i=0;i<numberNodes2_;i++) {
2365      int depth =  statistics_[i]->depth(); 
2366      int way =  statistics_[i]->way(); 
2367      double value = statistics_[i]->value(); 
2368      double startingObjective =  statistics_[i]->startingObjective(); 
2369      int startingInfeasibility = statistics_[i]->startingInfeasibility(); 
2370      double endingObjective = statistics_[i]->endingObjective(); 
2371      int endingInfeasibility = statistics_[i]->endingInfeasibility(); 
2372      maxDepth = CoinMax(depth,maxDepth);
2373      // Only for completed
2374      averageNumberIterations1 += statistics_[i]->numberIterations();
2375      averageValue += value;
2376      if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) {
2377        numberSolutions++;
2378        averageSolutionDepth += depth;
2379      }
2380      if (endingObjective==COIN_DBL_MAX) {
2381        numberCutoff++;
2382        averageCutoffDepth += depth;
2383        if (way<0) {
2384          numberDown++;
2385          numberCutoffDown++;
2386          if (way==-1)
2387            numberFirstDown++;
2388        } else {
2389          numberUp++;
2390          numberCutoffUp++;
2391          if (way==1)
2392            numberFirstUp++;
2393        }
2394      } else {
2395        averageSolvedDepth += depth;
2396        if (way<0) {
2397          numberDown++;
2398          averageInfDown += startingInfeasibility-endingInfeasibility;
2399          averageObjDown += endingObjective-startingObjective;
2400          if (way==-1)
2401            numberFirstDown++;
2402        } else {
2403          numberUp++;
2404          averageInfUp += startingInfeasibility-endingInfeasibility;
2405          averageObjUp += endingObjective-startingObjective;
2406          if (way==1)
2407            numberFirstUp++;
2408        }
2409      }
2410    }
2411    // Now print
2412    if (numberSolutions)
2413      averageSolutionDepth /= (double) numberSolutions;
2414    int numberSolved = numberNodes2_-numberCutoff;
2415    double averageNumberIterations2=numberIterations_-averageNumberIterations1
2416      -numberIterationsAtContinuous;
2417    if(numberCutoff) {
2418      averageCutoffDepth /= (double) numberCutoff;
2419      averageNumberIterations2 /= (double) numberCutoff;
2420    }
2421    if (numberNodes2_) 
2422      averageValue /= (double) numberNodes2_;
2423    if (numberSolved) {
2424      averageNumberIterations1 /= (double) numberSolved;
2425      averageSolvedDepth /= (double) numberSolved;
2426    }
2427    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
2428           numberSolutions,averageSolutionDepth);
2429    printf("average value of variable being branched on was %g\n",
2430           averageValue);
2431    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
2432           numberCutoff,averageCutoffDepth,averageNumberIterations2);
2433    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
2434           numberSolved,averageSolvedDepth,averageNumberIterations1);
2435    if (numberDown) {
2436      averageInfDown /= (double) numberDown;
2437      averageObjDown /= (double) numberDown;
2438    }
2439    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
2440           numberDown,numberFirstDown,numberDown-numberFirstDown,numberCutoffDown,
2441           averageInfDown,averageObjDown);
2442    if (numberUp) {
2443      averageInfUp /= (double) numberUp;
2444      averageObjUp /= (double) numberUp;
2445    }
2446    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
2447           numberUp,numberFirstUp,numberUp-numberFirstUp,numberCutoffUp,
2448           averageInfUp,averageObjUp);
2449    for ( i=0;i<numberNodes2_;i++) 
2450      delete statistics_[i];
2451    delete [] statistics_;
2452    statistics_=NULL;
2453    maximumStatistics_=0;
2454    delete [] lookup;
2455  }
2456/*
2457  If we think we have a solution, restore and confirm it with a call to
2458  setBestSolution().  We need to reset the cutoff value so as not to fathom
2459  the solution on bounds.  Note that calling setBestSolution( ..., true)
2460  leaves the continuousSolver_ bounds vectors fixed at the solution value.
2461
2462  Running resolve() here is a failsafe --- setBestSolution has already
2463  reoptimised using the continuousSolver_. If for some reason we fail to
2464  prove optimality, run the problem again after instructing the solver to
2465  tell us more.
2466
2467  If all looks good, replace solver_ with continuousSolver_, so that the
2468  outside world will be able to obtain information about the solution using
2469  public methods.
2470*/
2471  if (bestSolution_&&(solverCharacteristics_->solverType()<2||solverCharacteristics_->solverType()==4)) 
2472  { setCutoff(1.0e50) ; // As best solution should be worse than cutoff
2473    phase_=5;
2474    double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2475    if ((specialOptions_&4)!=0)
2476      bestObjective_ += 100.0*increment+1.0e-3; // only set if we are going to solve
2477    setBestSolution(CBC_END_SOLUTION,bestObjective_,bestSolution_,true) ;
2478    continuousSolver_->resolve() ;
2479    if (!continuousSolver_->isProvenOptimal())
2480    { continuousSolver_->messageHandler()->setLogLevel(2) ;
2481      continuousSolver_->initialSolve() ; }
2482    delete solver_ ;
2483    // above deletes solverCharacteristics_
2484    solverCharacteristics_ = NULL;
2485    solver_ = continuousSolver_ ;
2486    setPointers(solver_);
2487    continuousSolver_ = NULL ; }
2488/*
2489  Clean up dangling objects. continuousSolver_ may already be toast.
2490*/
2491  delete lastws ;
2492  delete [] whichGenerator_ ;
2493  whichGenerator_=NULL;
2494  delete [] lowerBefore ;
2495  delete [] upperBefore ;
2496  delete [] walkback_ ;
2497  walkback_ = NULL ;
2498  delete [] addedCuts_ ;
2499  addedCuts_ = NULL ;
2500  //delete persistentInfo;
2501  // Get rid of characteristics
2502  solverCharacteristics_=NULL;
2503  if (continuousSolver_)
2504  { delete continuousSolver_ ;
2505    continuousSolver_ = NULL ; }
2506/*
2507  Destroy global cuts by replacing with an empty OsiCuts object.
2508*/
2509  globalCuts_= OsiCuts() ;
2510  if (!bestSolution_) {
2511    // make sure lp solver is infeasible
2512    int numberColumns = solver_->getNumCols();
2513    const double * columnLower = solver_->getColLower();
2514    int iColumn;
2515    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2516      if (solver_->isInteger(iColumn)) 
2517        solver_->setColUpper(iColumn,columnLower[iColumn]);
2518    }
2519    solver_->initialSolve();
2520  }
2521  if (strategy_&&strategy_->preProcessState()>0) {
2522    // undo preprocessing
2523    CglPreProcess * process = strategy_->process();
2524    assert (process);
2525    int n = originalSolver->getNumCols();
2526    if (bestSolution_) {
2527      delete [] bestSolution_;
2528      bestSolution_ = new double [n];
2529      process->postProcess(*solver_);
2530    }
2531    strategy_->deletePreProcess();
2532    // Solution now back in originalSolver
2533    delete solver_;
2534    solver_=originalSolver;
2535    if (bestSolution_)
2536      memcpy(bestSolution_,solver_->getColSolution(),n*sizeof(double));
2537    // put back original objects if there were any
2538    if (originalObject) {
2539      int iColumn;
2540      assert (ownObjects_);
2541      for (iColumn=0;iColumn<numberObjects_;iColumn++) 
2542        delete object_[iColumn];
2543      delete [] object_;
2544      numberObjects_ = numberOriginalObjects;
2545      object_=originalObject;
2546      delete [] integerVariable_;
2547      numberIntegers_=0;
2548      for (iColumn=0;iColumn<n;iColumn++) {
2549        if (solver_->isInteger(iColumn))
2550          numberIntegers_++;
2551      }
2552      integerVariable_ = new int[numberIntegers_];
2553      numberIntegers_=0;
2554      for (iColumn=0;iColumn<n;iColumn++) {
2555        if (solver_->isInteger(iColumn))
2556            integerVariable_[numberIntegers_++]=iColumn;
2557      }
2558    }
2559  }
2560#ifdef CLP_QUICK_OPTIONS
2561  {
2562    OsiClpSolverInterface * clpSolver
2563      = dynamic_cast<OsiClpSolverInterface *> (solver_);
2564    if (clpSolver) {
2565      // Try and re-use regions   
2566      ClpSimplex * simplex = clpSolver->getModelPtr();
2567      simplex->setPersistenceFlag(0);
2568      clpSolver->deleteScaleFactors();
2569      clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~131072));
2570      simplex->setSpecialOptions(simplex->specialOptions()&(~131072));
2571      simplex->setAlphaAccuracy(-1.0);
2572      //clpSolver->setSpecialOptions((clpSolver->specialOptions()&~128)|65536);
2573    }
2574  }
2575#endif
2576  return ;
2577 }
2578
2579
2580// Solve the initial LP relaxation
2581void 
2582CbcModel::initialSolve() 
2583{
2584  assert (solver_);
2585  assert (!solverCharacteristics_);
2586  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
2587  if (solverCharacteristics) {
2588    solverCharacteristics_ = solverCharacteristics;
2589  } else {
2590    // replace in solver
2591    OsiBabSolver defaultC;
2592    solver_->setAuxiliaryInfo(&defaultC);
2593    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
2594  }
2595  solverCharacteristics_->setSolver(solver_);
2596  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2597  solver_->initialSolve();
2598  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
2599  // But set up so Jon Lee will be happy
2600  status_=-1;
2601  secondaryStatus_ = -1;
2602  originalContinuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
2603  delete [] continuousSolution_;
2604  continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
2605                                             solver_->getNumCols());
2606  setPointers(solver_);
2607  solverCharacteristics_ = NULL;
2608}
2609
2610/*! \brief Get an empty basis object
2611
2612  Return an empty CoinWarmStartBasis object with the requested capacity,
2613  appropriate for the current solver. The object is cloned from the object
2614  cached as emptyWarmStart_. If there is no cached object, the routine
2615  queries the solver for a warm start object, empties it, and caches the
2616  result.
2617*/
2618
2619CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
2620
2621{ CoinWarmStartBasis *emptyBasis ;
2622/*
2623  Acquire an empty basis object, if we don't yet have one.
2624*/
2625  if (emptyWarmStart_ == 0)
2626  { if (solver_ == 0)
2627    { throw CoinError("Cannot construct basis without solver!",
2628                      "getEmptyBasis","CbcModel") ; }
2629    emptyBasis =
2630        dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
2631    if (emptyBasis == 0)
2632    { throw CoinError(
2633        "Solver does not appear to use a basis-oriented warm start.",
2634        "getEmptyBasis","CbcModel") ; }
2635    emptyBasis->setSize(0,0) ;
2636    emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ; }
2637/*
2638  Clone the empty basis object, resize it as requested, and return.
2639*/
2640  emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
2641  assert(emptyBasis) ;
2642  if (ns != 0 || na != 0) emptyBasis->setSize(ns,na) ;
2643
2644  return (emptyBasis) ; }
2645   
2646
2647/** Default Constructor
2648
2649  Creates an empty model without an associated solver.
2650*/
2651CbcModel::CbcModel() 
2652
2653:
2654  solver_(NULL),
2655  ourSolver_(true),
2656  continuousSolver_(NULL),
2657  referenceSolver_(NULL),
2658  defaultHandler_(true),
2659  emptyWarmStart_(NULL),
2660  bestObjective_(COIN_DBL_MAX),
2661  bestPossibleObjective_(COIN_DBL_MAX),
2662  sumChangeObjective1_(0.0),
2663  sumChangeObjective2_(0.0),
2664  bestSolution_(NULL),
2665  currentSolution_(NULL),
2666  testSolution_(NULL),
2667  minimumDrop_(1.0e-4),
2668  numberSolutions_(0),
2669  stateOfSearch_(0),
2670  hotstartSolution_(NULL),
2671  hotstartPriorities_(NULL),
2672  numberHeuristicSolutions_(0),
2673  numberNodes_(0),
2674  numberNodes2_(0),
2675  numberIterations_(0),
2676  status_(-1),
2677  secondaryStatus_(-1),
2678  numberIntegers_(0),
2679  numberRowsAtContinuous_(0),
2680  maximumNumberCuts_(0),
2681  phase_(0),
2682  currentNumberCuts_(0),
2683  maximumDepth_(0),
2684  walkback_(NULL),
2685  addedCuts_(NULL),
2686  nextRowCut_(NULL),
2687  currentNode_(NULL),
2688  integerVariable_(NULL),
2689  integerInfo_(NULL),
2690  continuousSolution_(NULL),
2691  usedInSolution_(NULL),
2692  specialOptions_(0),
2693  subTreeModel_(NULL),
2694  numberStoppedSubTrees_(0),
2695  mutex_(NULL),
2696  presolve_(0),
2697  numberStrong_(5),
2698  numberBeforeTrust_(0),
2699  numberPenalties_(20),
2700  penaltyScaleFactor_(3.0),
2701  numberAnalyzeIterations_(0),
2702  analyzeResults_(NULL),
2703  numberInfeasibleNodes_(0),
2704  problemType_(0),
2705  printFrequency_(0),
2706  numberCutGenerators_(0),
2707  generator_(NULL),
2708  virginGenerator_(NULL),
2709  numberHeuristics_(0),
2710  heuristic_(NULL),
2711  lastHeuristic_(NULL),
2712  eventHandler_(0),
2713  numberObjects_(0),
2714  object_(NULL),
2715  ownObjects_(true),
2716  originalColumns_(NULL),
2717  howOftenGlobalScan_(1),
2718  numberGlobalViolations_(0),
2719  continuousObjective_(COIN_DBL_MAX),
2720  originalContinuousObjective_(COIN_DBL_MAX),
2721  continuousInfeasibilities_(INT_MAX),
2722  maximumCutPassesAtRoot_(20),
2723  maximumCutPasses_(10),
2724  preferredWay_(0),
2725  currentPassNumber_(0),
2726  maximumWhich_(1000),
2727  whichGenerator_(NULL),
2728  maximumStatistics_(0),
2729  statistics_(NULL),
2730  maximumDepthActual_(0),
2731  numberDJFixed_(0.0),
2732  probingInfo_(NULL),
2733  numberFixedAtRoot_(0),
2734  numberFixedNow_(0),
2735  stoppedOnGap_(false),
2736  eventHappened_(false),
2737  numberLongStrong_(0),
2738  numberOldActiveCuts_(0),
2739  numberNewCuts_(0),
2740  sizeMiniTree_(0),
2741  searchStrategy_(-1),
2742  numberStrongIterations_(0),
2743  resolveAfterTakeOffCuts_(true),
2744  numberThreads_(0),
2745  threadMode_(0)
2746{
2747  memset(intParam_,0,sizeof(intParam_));
2748  intParam_[CbcMaxNumNode] = 2147483647;
2749  intParam_[CbcMaxNumSol] = 9999999;
2750  intParam_[CbcFathomDiscipline] = 0;
2751
2752  dblParam_[CbcIntegerTolerance] = 1e-6;
2753  dblParam_[CbcInfeasibilityWeight] = 0.0;
2754  dblParam_[CbcCutoffIncrement] = 1e-5;
2755  dblParam_[CbcAllowableGap] = 1.0e-10;
2756  dblParam_[CbcAllowableFractionGap] = 0.0;
2757  dblParam_[CbcMaximumSeconds] = 1.0e100;
2758  dblParam_[CbcCurrentCutoff] = 1.0e100;
2759  dblParam_[CbcOptimizationDirection] = 1.0;
2760  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
2761  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
2762  dblParam_[CbcStartSeconds] = 0.0;
2763  strongInfo_[0]=0;
2764  strongInfo_[1]=0;
2765  strongInfo_[2]=0;
2766  solverCharacteristics_ = NULL;
2767  nodeCompare_=new CbcCompareDefault();;
2768  problemFeasibility_=new CbcFeasibilityBase();
2769  tree_= new CbcTree();
2770  branchingMethod_=NULL;
2771  cutModifier_=NULL;
2772  strategy_=NULL;
2773  parentModel_=NULL;
2774  cbcColLower_ = NULL;
2775  cbcColUpper_ = NULL;
2776  cbcRowLower_ = NULL;
2777  cbcRowUpper_ = NULL;
2778  cbcColSolution_ = NULL;
2779  cbcRowPrice_ = NULL;
2780  cbcReducedCost_ = NULL;
2781  cbcRowActivity_ = NULL;
2782  appData_=NULL;
2783  handler_ = new CoinMessageHandler();
2784  handler_->setLogLevel(2);
2785  messages_ = CbcMessage();
2786  eventHandler_ = new CbcEventHandler() ;
2787}
2788
2789/** Constructor from solver.
2790
2791  Creates a model complete with a clone of the solver passed as a parameter.
2792*/
2793
2794CbcModel::CbcModel(const OsiSolverInterface &rhs)
2795:
2796  continuousSolver_(NULL),
2797  referenceSolver_(NULL),
2798  defaultHandler_(true),
2799  emptyWarmStart_(NULL),
2800  bestObjective_(COIN_DBL_MAX),
2801  bestPossibleObjective_(COIN_DBL_MAX),
2802  sumChangeObjective1_(0.0),
2803  sumChangeObjective2_(0.0),
2804  minimumDrop_(1.0e-4),
2805  numberSolutions_(0),
2806  stateOfSearch_(0),
2807  hotstartSolution_(NULL),
2808  hotstartPriorities_(NULL),
2809  numberHeuristicSolutions_(0),
2810  numberNodes_(0),
2811  numberNodes2_(0),
2812  numberIterations_(0),
2813  status_(-1),
2814  secondaryStatus_(-1),
2815  numberRowsAtContinuous_(0),
2816  maximumNumberCuts_(0),
2817  phase_(0),
2818  currentNumberCuts_(0),
2819  maximumDepth_(0),
2820  walkback_(NULL),
2821  addedCuts_(NULL),
2822  nextRowCut_(NULL),
2823  currentNode_(NULL),
2824  integerInfo_(NULL),
2825  specialOptions_(0),
2826  subTreeModel_(NULL),
2827  numberStoppedSubTrees_(0),
2828  mutex_(NULL),
2829  presolve_(0),
2830  numberStrong_(5),
2831  numberBeforeTrust_(0),
2832  numberPenalties_(20),
2833  penaltyScaleFactor_(3.0),
2834  numberAnalyzeIterations_(0),
2835  analyzeResults_(NULL),
2836  numberInfeasibleNodes_(0),
2837  problemType_(0),
2838  printFrequency_(0),
2839  numberCutGenerators_(0),
2840  generator_(NULL),
2841  virginGenerator_(NULL),
2842  numberHeuristics_(0),
2843  heuristic_(NULL),
2844  lastHeuristic_(NULL),
2845  eventHandler_(0),
2846  numberObjects_(0),
2847  object_(NULL),
2848  ownObjects_(true),
2849  originalColumns_(NULL),
2850  howOftenGlobalScan_(1),
2851  numberGlobalViolations_(0),
2852  continuousObjective_(COIN_DBL_MAX),
2853  originalContinuousObjective_(COIN_DBL_MAX),
2854  continuousInfeasibilities_(INT_MAX),
2855  maximumCutPassesAtRoot_(20),
2856  maximumCutPasses_(10),
2857  preferredWay_(0),
2858  currentPassNumber_(0),
2859  maximumWhich_(1000),
2860  whichGenerator_(NULL),
2861  maximumStatistics_(0),
2862  statistics_(NULL),
2863  maximumDepthActual_(0),
2864  numberDJFixed_(0.0),
2865  probingInfo_(NULL),
2866  numberFixedAtRoot_(0),
2867  numberFixedNow_(0),
2868  stoppedOnGap_(false),
2869  eventHappened_(false),
2870  numberLongStrong_(0),
2871  numberOldActiveCuts_(0),
2872  numberNewCuts_(0),
2873  sizeMiniTree_(0),
2874  searchStrategy_(-1),
2875  numberStrongIterations_(0),
2876  resolveAfterTakeOffCuts_(true),
2877  numberThreads_(0),
2878  threadMode_(0)
2879{
2880  memset(intParam_,0,sizeof(intParam_));
2881  intParam_[CbcMaxNumNode] = 2147483647;
2882  intParam_[CbcMaxNumSol] = 9999999;
2883  intParam_[CbcFathomDiscipline] = 0;
2884
2885  dblParam_[CbcIntegerTolerance] = 1e-6;
2886  dblParam_[CbcInfeasibilityWeight] = 0.0;
2887  dblParam_[CbcCutoffIncrement] = 1e-5;
2888  dblParam_[CbcAllowableGap] = 1.0e-10;
2889  dblParam_[CbcAllowableFractionGap] = 0.0;
2890  dblParam_[CbcMaximumSeconds] = 1.0e100;
2891  dblParam_[CbcCurrentCutoff] = 1.0e100;
2892  dblParam_[CbcOptimizationDirection] = 1.0;
2893  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
2894  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
2895  dblParam_[CbcStartSeconds] = 0.0;
2896  strongInfo_[0]=0;
2897  strongInfo_[1]=0;
2898  strongInfo_[2]=0;
2899  solverCharacteristics_ = NULL;
2900
2901  nodeCompare_=new CbcCompareDefault();;
2902  problemFeasibility_=new CbcFeasibilityBase();
2903  tree_= new CbcTree();
2904  branchingMethod_=NULL;
2905  cutModifier_=NULL;
2906  strategy_=NULL;
2907  parentModel_=NULL;
2908  appData_=NULL;
2909  handler_ = new CoinMessageHandler();
2910  handler_->setLogLevel(2);
2911  messages_ = CbcMessage();
2912  eventHandler_ = new CbcEventHandler() ;
2913  solver_ = rhs.clone();
2914  referenceSolver_ = solver_->clone();
2915  ourSolver_ = true ;
2916  cbcColLower_ = NULL;
2917  cbcColUpper_ = NULL;
2918  cbcRowLower_ = NULL;
2919  cbcRowUpper_ = NULL;
2920  cbcColSolution_ = NULL;
2921  cbcRowPrice_ = NULL;
2922  cbcReducedCost_ = NULL;
2923  cbcRowActivity_ = NULL;
2924
2925  // Initialize solution and integer variable vectors
2926  bestSolution_ = NULL; // to say no solution found
2927  numberIntegers_=0;
2928  int numberColumns = solver_->getNumCols();
2929  int iColumn;
2930  if (numberColumns) {
2931    // Space for current solution
2932    currentSolution_ = new double[numberColumns];
2933    continuousSolution_ = new double[numberColumns];
2934    usedInSolution_ = new int[numberColumns];
2935    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2936      if( solver_->isInteger(iColumn)) 
2937        numberIntegers_++;
2938    }
2939  } else {
2940    // empty model
2941    currentSolution_=NULL;
2942    continuousSolution_=NULL;
2943    usedInSolution_=NULL;
2944  }
2945  testSolution_=currentSolution_;
2946  if (numberIntegers_) {
2947    integerVariable_ = new int [numberIntegers_];
2948    numberIntegers_=0;
2949    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2950      if( solver_->isInteger(iColumn)) 
2951        integerVariable_[numberIntegers_++]=iColumn;
2952    }
2953  } else {
2954    integerVariable_ = NULL;
2955  }
2956}
2957
2958/*
2959  Assign a solver to the model (model assumes ownership)
2960
2961  The integer variable vector is initialized if it's not already present.
2962  If deleteSolver then current solver deleted (if model owned)
2963
2964  Assuming ownership matches usage in OsiSolverInterface
2965  (cf. assignProblem, loadProblem).
2966
2967  TODO: What to do about solver parameters? A simple copy likely won't do it,
2968        because the SI must push the settings into the underlying solver. In
2969        the context of switching solvers in cbc, this means that command line
2970        settings will get lost. Stash the command line somewhere and reread it
2971        here, maybe?
2972 
2973  TODO: More generally, how much state should be transferred from the old
2974        solver to the new solver? Best perhaps to see how usage develops.
2975        What's done here mimics the CbcModel(OsiSolverInterface) constructor.
2976*/
2977void
2978CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver)
2979
2980{
2981  // resize best solution if exists
2982  if (bestSolution_) {
2983    int nOld = solver_->getNumCols();
2984    int nNew = solver->getNumCols();
2985    if (nNew>nOld) {
2986      double * temp = new double[nNew];
2987      memcpy(temp,bestSolution_,nOld*sizeof(double));
2988      memset(temp+nOld,0,(nNew-nOld)*sizeof(double));
2989      delete [] bestSolution_;
2990      bestSolution_=temp;
2991    }
2992  }
2993  // Keep the current message level for solver (if solver exists)
2994  if (solver_)
2995    solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
2996
2997  if (ourSolver_&&deleteSolver) delete solver_ ;
2998  solver_ = solver;
2999  solver = NULL ;
3000  ourSolver_ = true ;
3001/*
3002  Basis information is solver-specific.
3003*/
3004  if (emptyWarmStart_)
3005  { delete emptyWarmStart_  ;
3006    emptyWarmStart_ = 0 ; }
3007/*
3008  Initialize integer variable vector.
3009*/
3010  numberIntegers_=0;
3011  int numberColumns = solver_->getNumCols();
3012  int iColumn;
3013  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3014    if( solver_->isInteger(iColumn)) 
3015      numberIntegers_++;
3016  }
3017  if (numberIntegers_) {
3018    delete [] integerVariable_;
3019    integerVariable_ = new int [numberIntegers_];
3020    numberIntegers_=0;
3021    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3022      if( solver_->isInteger(iColumn)) 
3023        integerVariable_[numberIntegers_++]=iColumn;
3024    }
3025  } else {
3026    integerVariable_ = NULL;
3027  }
3028
3029  return ;
3030}
3031
3032// Copy constructor.
3033
3034CbcModel::CbcModel(const CbcModel & rhs, bool noTree)
3035:
3036  continuousSolver_(NULL),
3037  referenceSolver_(NULL),
3038  defaultHandler_(rhs.defaultHandler_),
3039  emptyWarmStart_(NULL),
3040  bestObjective_(rhs.bestObjective_),
3041  bestPossibleObjective_(rhs.bestPossibleObjective_),
3042  sumChangeObjective1_(rhs.sumChangeObjective1_),
3043  sumChangeObjective2_(rhs.sumChangeObjective2_),
3044  minimumDrop_(rhs.minimumDrop_),
3045  numberSolutions_(rhs.numberSolutions_),
3046  stateOfSearch_(rhs.stateOfSearch_),
3047  numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
3048  numberNodes_(rhs.numberNodes_),
3049  numberNodes2_(rhs.numberNodes2_),
3050  numberIterations_(rhs.numberIterations_),
3051  status_(rhs.status_),
3052  secondaryStatus_(rhs.secondaryStatus_),
3053  specialOptions_(rhs.specialOptions_),
3054  subTreeModel_(rhs.subTreeModel_),
3055  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
3056  mutex_(NULL),
3057  presolve_(rhs.presolve_),
3058  numberStrong_(rhs.numberStrong_),
3059  numberBeforeTrust_(rhs.numberBeforeTrust_),
3060  numberPenalties_(rhs.numberPenalties_),
3061  penaltyScaleFactor_(rhs.penaltyScaleFactor_),
3062  numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
3063  analyzeResults_(NULL),
3064  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
3065  problemType_(rhs.problemType_),
3066  printFrequency_(rhs.printFrequency_),
3067  howOftenGlobalScan_(rhs.howOftenGlobalScan_),
3068  numberGlobalViolations_(rhs.numberGlobalViolations_),
3069  continuousObjective_(rhs.continuousObjective_),
3070  originalContinuousObjective_(rhs.originalContinuousObjective_),
3071  continuousInfeasibilities_(rhs.continuousInfeasibilities_),
3072  maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
3073  maximumCutPasses_( rhs.maximumCutPasses_),
3074  preferredWay_(rhs.preferredWay_),
3075  currentPassNumber_(rhs.currentPassNumber_),
3076  maximumWhich_(rhs.maximumWhich_),
3077  whichGenerator_(NULL),
3078  maximumStatistics_(0),
3079  statistics_(NULL),
3080  maximumDepthActual_(0),
3081  numberDJFixed_(0.0),
3082  probingInfo_(NULL),
3083  numberFixedAtRoot_(rhs.numberFixedAtRoot_),
3084  numberFixedNow_(rhs.numberFixedNow_),
3085  stoppedOnGap_(rhs.stoppedOnGap_),
3086  eventHappened_(rhs.eventHappened_),
3087  numberLongStrong_(rhs.numberLongStrong_),
3088  numberOldActiveCuts_(rhs.numberOldActiveCuts_),
3089  numberNewCuts_(rhs.numberNewCuts_),
3090  sizeMiniTree_(rhs.sizeMiniTree_),
3091  searchStrategy_(rhs.searchStrategy_),
3092  numberStrongIterations_(rhs.numberStrongIterations_),
3093  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
3094  numberThreads_(rhs.numberThreads_),
3095  threadMode_(rhs.threadMode_)
3096{
3097  memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
3098  memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
3099  strongInfo_[0]=rhs.strongInfo_[0];
3100  strongInfo_[1]=rhs.strongInfo_[1];
3101  strongInfo_[2]=rhs.strongInfo_[2];
3102  solverCharacteristics_ = NULL;
3103  if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
3104  if (defaultHandler_) {
3105    handler_ = new CoinMessageHandler();
3106    handler_->setLogLevel(2);
3107  } else {
3108    handler_ = rhs.handler_;
3109  }
3110  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3111  numberCutGenerators_ = rhs.numberCutGenerators_;
3112  if (numberCutGenerators_) {
3113    generator_ = new CbcCutGenerator * [numberCutGenerators_];
3114    virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
3115    int i;
3116    for (i=0;i<numberCutGenerators_;i++) {
3117      generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
3118      virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
3119    }
3120  } else {
3121    generator_=NULL;
3122    virginGenerator_=NULL;
3123  }
3124  if (!noTree)
3125    globalCuts_ = rhs.globalCuts_;
3126  numberHeuristics_ = rhs.numberHeuristics_;
3127  if (numberHeuristics_) {
3128    heuristic_ = new CbcHeuristic * [numberHeuristics_];
3129    int i;
3130    for (i=0;i<numberHeuristics_;i++) {
3131      heuristic_[i]=rhs.heuristic_[i]->clone();
3132    }
3133  } else {
3134    heuristic_=NULL;
3135  }
3136  lastHeuristic_ = NULL;
3137  if (rhs.eventHandler_)
3138  { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; }
3139  else
3140  { eventHandler_ = NULL ; }
3141  ownObjects_ = rhs.ownObjects_;
3142  if (ownObjects_) {
3143    numberObjects_=rhs.numberObjects_;
3144    if (numberObjects_) {
3145      object_ = new OsiObject * [numberObjects_];
3146      int i;
3147      for (i=0;i<numberObjects_;i++) {
3148        object_[i]=(rhs.object_[i])->clone();
3149        CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
3150        assert (obj);
3151        obj->setModel(this);
3152      }
3153    } else {
3154      object_=NULL;
3155    }
3156  } else {
3157    // assume will be redone
3158    numberObjects_=0;
3159    object_=NULL;
3160  }
3161  if (rhs.referenceSolver_)
3162    referenceSolver_ = rhs.referenceSolver_->clone();
3163  else
3164    referenceSolver_=NULL;
3165  if (!noTree||!rhs.continuousSolver_)
3166    solver_ = rhs.solver_->clone();
3167  else
3168    solver_ = rhs.continuousSolver_->clone();
3169  if (rhs.originalColumns_) {
3170    int numberColumns = solver_->getNumCols();
3171    originalColumns_= new int [numberColumns];
3172    memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
3173  } else {
3174    originalColumns_=NULL;
3175  }
3176  if (maximumWhich_&&rhs.whichGenerator_)
3177    whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
3178  nodeCompare_=rhs.nodeCompare_->clone();
3179  problemFeasibility_=rhs.problemFeasibility_->clone();
3180  tree_= rhs.tree_->clone();
3181  if (rhs.branchingMethod_)
3182    branchingMethod_=rhs.branchingMethod_->clone();
3183  else
3184    branchingMethod_=NULL;
3185  if (rhs.cutModifier_)
3186    cutModifier_=rhs.cutModifier_->clone();
3187  else
3188    cutModifier_=NULL;
3189  cbcColLower_ = NULL;
3190  cbcColUpper_ = NULL;
3191  cbcRowLower_ = NULL;
3192  cbcRowUpper_ = NULL;
3193  cbcColSolution_ = NULL;
3194  cbcRowPrice_ = NULL;
3195  cbcReducedCost_ = NULL;
3196  cbcRowActivity_ = NULL;
3197  if (rhs.strategy_)
3198    strategy_=rhs.strategy_->clone();
3199  else
3200    strategy_=NULL;
3201  parentModel_=rhs.parentModel_;
3202  appData_=rhs.appData_;
3203  messages_ = rhs.messages_;
3204  ourSolver_ = true ;
3205  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3206  numberIntegers_=rhs.numberIntegers_;
3207  if (numberIntegers_) {
3208    integerVariable_ = new int [numberIntegers_];
3209    memcpy(integerVariable_,rhs.integerVariable_,numberIntegers_*sizeof(int));
3210    integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
3211  } else {
3212    integerVariable_ = NULL;
3213    integerInfo_=NULL;
3214  }
3215  if (rhs.hotstartSolution_) {
3216    int numberColumns = solver_->getNumCols();
3217    hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
3218    hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
3219  } else {
3220    hotstartSolution_ = NULL;
3221    hotstartPriorities_ =NULL;
3222  }
3223  if (rhs.bestSolution_&&!noTree) {
3224    int numberColumns = solver_->getNumCols();
3225    bestSolution_ = new double[numberColumns];
3226    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
3227  } else {
3228    bestSolution_=NULL;
3229  }
3230  if (!noTree) {
3231    int numberColumns = solver_->getNumCols();
3232    currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns);
3233    continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns);
3234    usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns);
3235  } else {
3236    currentSolution_=NULL;
3237    continuousSolution_=NULL;
3238    usedInSolution_=NULL;
3239  }
3240  testSolution_=currentSolution_;
3241  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
3242  maximumNumberCuts_=rhs.maximumNumberCuts_;
3243  phase_ = rhs.phase_;
3244  currentNumberCuts_=rhs.currentNumberCuts_;
3245  maximumDepth_= rhs.maximumDepth_;
3246  if (noTree) {
3247    bestObjective_ = COIN_DBL_MAX;
3248    numberSolutions_ =0;
3249    stateOfSearch_= 0;
3250    numberHeuristicSolutions_=0;
3251    numberNodes_=0;
3252    numberNodes2_=0;
3253    numberIterations_=0;
3254    status_=0;
3255    subTreeModel_=NULL;
3256    numberStoppedSubTrees_=0;
3257    continuousObjective_=COIN_DBL_MAX;
3258    originalContinuousObjective_=COIN_DBL_MAX;
3259    continuousInfeasibilities_=INT_MAX;
3260    maximumNumberCuts_=0;
3261    tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_);
3262    bestPossibleObjective_ = COIN_DBL_MAX;
3263  }
3264  // These are only used as temporary arrays so need not be filled
3265  if (maximumNumberCuts_) {
3266    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3267  } else {
3268    addedCuts_ = NULL;
3269  }
3270  nextRowCut_ = NULL;
3271  currentNode_ = NULL;
3272  if (maximumDepth_)
3273    walkback_ = new CbcNodeInfo * [maximumDepth_];
3274  else
3275    walkback_ = NULL;
3276  synchronizeModel();
3277}
3278 
3279// Assignment operator
3280CbcModel & 
3281CbcModel::operator=(const CbcModel& rhs)
3282{
3283  if (this!=&rhs) {
3284    gutsOfDestructor();
3285    if (defaultHandler_)
3286    { delete handler_;
3287      handler_ = NULL; }
3288    defaultHandler_ = rhs.defaultHandler_;
3289    if (defaultHandler_)
3290    { handler_ = new CoinMessageHandler();
3291      handler_->setLogLevel(2); }
3292    else
3293    { handler_ = rhs.handler_; }
3294    messages_ = rhs.messages_;
3295    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3296
3297    delete solver_;
3298    if (rhs.solver_)
3299    { solver_ = rhs.solver_->clone() ; }
3300    else
3301    { solver_ = 0 ; }
3302    ourSolver_ = true ;
3303    delete continuousSolver_ ;
3304    if (rhs.continuousSolver_)
3305    { solver_ = rhs.continuousSolver_->clone() ; }
3306    else
3307    { continuousSolver_ = 0 ; }
3308
3309    if (rhs.referenceSolver_)
3310    { solver_ = rhs.referenceSolver_->clone() ; }
3311    else
3312    { referenceSolver_ = NULL ; }
3313
3314    delete emptyWarmStart_ ;
3315    if (rhs.emptyWarmStart_)
3316    { emptyWarmStart_ = rhs.emptyWarmStart_->clone() ; }
3317    else
3318    { emptyWarmStart_ = 0 ; }
3319
3320    bestObjective_ = rhs.bestObjective_;
3321    bestPossibleObjective_=rhs.bestPossibleObjective_;
3322    sumChangeObjective1_=rhs.sumChangeObjective1_;
3323    sumChangeObjective2_=rhs.sumChangeObjective2_;
3324    delete [] bestSolution_;
3325    if (rhs.bestSolution_) {
3326      int numberColumns = rhs.getNumCols();
3327      bestSolution_ = new double[numberColumns];
3328      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
3329    } else {
3330      bestSolution_=NULL;
3331    }
3332    int numberColumns = solver_->getNumCols();
3333    currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns);
3334    continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns);
3335    usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns);
3336    testSolution_=currentSolution_;
3337    minimumDrop_ = rhs.minimumDrop_;
3338    numberSolutions_=rhs.numberSolutions_;
3339    stateOfSearch_= rhs.stateOfSearch_;
3340    numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
3341    numberNodes_ = rhs.numberNodes_;
3342    numberNodes2_ = rhs.numberNodes2_;
3343    numberIterations_ = rhs.numberIterations_;
3344    status_ = rhs.status_;
3345    secondaryStatus_ = rhs.secondaryStatus_;
3346    specialOptions_ = rhs.specialOptions_;
3347    subTreeModel_ = rhs.subTreeModel_;
3348    numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
3349    mutex_ = NULL;
3350    presolve_ = rhs.presolve_;
3351    numberStrong_ = rhs.numberStrong_;
3352    numberBeforeTrust_ = rhs.numberBeforeTrust_;
3353    numberPenalties_ = rhs.numberPenalties_;
3354    penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
3355    numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
3356    delete [] analyzeResults_;
3357    analyzeResults_ = NULL;
3358    numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
3359    problemType_ = rhs.problemType_;
3360    printFrequency_ = rhs.printFrequency_;
3361    howOftenGlobalScan_=rhs.howOftenGlobalScan_;
3362    numberGlobalViolations_=rhs.numberGlobalViolations_;
3363    continuousObjective_=rhs.continuousObjective_;
3364    originalContinuousObjective_ = rhs.originalContinuousObjective_;
3365    continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
3366    maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
3367    maximumCutPasses_ = rhs.maximumCutPasses_;
3368    preferredWay_ = rhs.preferredWay_;
3369    currentPassNumber_ = rhs.currentPassNumber_;
3370    memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
3371    memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
3372    globalCuts_ = rhs.globalCuts_;
3373    int i;
3374    for (i=0;i<numberCutGenerators_;i++) {
3375      delete generator_[i];
3376      delete virginGenerator_[i];
3377    }
3378    delete [] generator_;
3379    delete [] virginGenerator_;
3380    delete [] heuristic_;
3381    maximumWhich_ = rhs.maximumWhich_;
3382    delete [] whichGenerator_;
3383    whichGenerator_ = NULL;
3384    if (maximumWhich_&&rhs.whichGenerator_)
3385      whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
3386    for (i=0;i<maximumStatistics_;i++)
3387      delete statistics_[i];
3388    delete [] statistics_;
3389    maximumStatistics_ = 0;
3390    statistics_ = NULL;
3391    delete probingInfo_;
3392    probingInfo_=NULL;
3393    numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
3394    numberFixedNow_ = rhs.numberFixedNow_;
3395    stoppedOnGap_ = rhs.stoppedOnGap_;
3396    eventHappened_ = rhs.eventHappened_;
3397    numberLongStrong_ = rhs.numberLongStrong_;
3398    numberOldActiveCuts_ = rhs.numberOldActiveCuts_;
3399    numberNewCuts_ = rhs.numberNewCuts_;
3400    resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_;
3401    numberThreads_ = rhs.numberThreads_;
3402    threadMode_ = rhs.threadMode_;
3403    sizeMiniTree_ = rhs.sizeMiniTree_;
3404    searchStrategy_ = rhs.searchStrategy_;
3405    numberStrongIterations_ = rhs.numberStrongIterations_;
3406    strongInfo_[0]=rhs.strongInfo_[0];
3407    strongInfo_[1]=rhs.strongInfo_[1];
3408    strongInfo_[2]=rhs.strongInfo_[2];
3409    solverCharacteristics_ = NULL;
3410    lastHeuristic_ = NULL;
3411    numberCutGenerators_ = rhs.numberCutGenerators_;
3412    if (numberCutGenerators_) {
3413      generator_ = new CbcCutGenerator * [numberCutGenerators_];
3414      virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
3415      int i;
3416      for (i=0;i<numberCutGenerators_;i++) {
3417        generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
3418        virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
3419      }
3420    } else {
3421      generator_=NULL;
3422      virginGenerator_=NULL;
3423    }
3424    numberHeuristics_ = rhs.numberHeuristics_;
3425    if (numberHeuristics_) {
3426      heuristic_ = new CbcHeuristic * [numberHeuristics_];
3427      memcpy(heuristic_,rhs.heuristic_,
3428             numberHeuristics_*sizeof(CbcHeuristic *));
3429    } else {
3430      heuristic_=NULL;
3431    }
3432    lastHeuristic_ = NULL;
3433    if (eventHandler_)
3434      delete eventHandler_ ;
3435    if (rhs.eventHandler_)
3436    { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; }
3437    else
3438    { eventHandler_ = NULL ; }
3439    if (ownObjects_) {
3440      for (i=0;i<numberObjects_;i++)
3441        delete object_[i];
3442      delete [] object_;
3443      numberObjects_=rhs.numberObjects_;
3444      if (numberObjects_) {
3445        object_ = new OsiObject * [numberObjects_];
3446        int i;
3447        for (i=0;i<numberObjects_;i++) 
3448          object_[i]=(rhs.object_[i])->clone();
3449      } else {
3450        object_=NULL;
3451    }
3452    } else {
3453      // assume will be redone
3454      numberObjects_=0;
3455      object_=NULL;
3456    }
3457    delete [] originalColumns_;
3458    if (rhs.originalColumns_) {
3459      int numberColumns = rhs.getNumCols();
3460      originalColumns_= new int [numberColumns];
3461      memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
3462    } else {
3463      originalColumns_=NULL;
3464    }
3465    nodeCompare_=rhs.nodeCompare_->clone();
3466    problemFeasibility_=rhs.problemFeasibility_->clone();
3467    delete tree_;
3468    tree_= rhs.tree_->clone();
3469    if (rhs.branchingMethod_)
3470      branchingMethod_=rhs.branchingMethod_->clone();
3471    else
3472      branchingMethod_=NULL;
3473    if (rhs.cutModifier_)
3474      cutModifier_=rhs.cutModifier_->clone();
3475    else
3476      cutModifier_=NULL;
3477    delete strategy_;
3478    if (rhs.strategy_)
3479      strategy_=rhs.strategy_->clone();
3480    else
3481      strategy_=NULL;
3482    parentModel_=rhs.parentModel_;
3483    appData_=rhs.appData_;
3484
3485    delete [] integerVariable_;
3486    numberIntegers_=rhs.numberIntegers_;
3487    if (numberIntegers_) {
3488      integerVariable_ = new int [numberIntegers_];
3489      memcpy(integerVariable_,rhs.integerVariable_,
3490             numberIntegers_*sizeof(int));
3491      integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
3492    } else {
3493      integerVariable_ = NULL;
3494      integerInfo_=NULL;
3495    }
3496    if (rhs.hotstartSolution_) {
3497      int numberColumns = solver_->getNumCols();
3498      hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
3499      hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
3500    } else {
3501      hotstartSolution_ = NULL;
3502      hotstartPriorities_ =NULL;
3503    }
3504    numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
3505    maximumNumberCuts_=rhs.maximumNumberCuts_;
3506    phase_ = rhs.phase_;
3507    currentNumberCuts_=rhs.currentNumberCuts_;
3508    maximumDepth_= rhs.maximumDepth_;
3509    delete [] addedCuts_;
3510    delete [] walkback_;
3511    // These are only used as temporary arrays so need not be filled
3512    if (maximumNumberCuts_) {
3513      addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3514    } else {
3515      addedCuts_ = NULL;
3516    }
3517    nextRowCut_ = NULL;
3518    currentNode_ = NULL;
3519    if (maximumDepth_)
3520      walkback_ = new CbcNodeInfo * [maximumDepth_];
3521    else
3522      walkback_ = NULL;
3523    synchronizeModel();
3524    cbcColLower_ = NULL;
3525    cbcColUpper_ = NULL;
3526    cbcRowLower_ = NULL;
3527    cbcRowUpper_ = NULL;
3528    cbcColSolution_ = NULL;
3529    cbcRowPrice_ = NULL;
3530    cbcReducedCost_ = NULL;
3531    cbcRowActivity_ = NULL;
3532  }
3533  return *this;
3534}
3535 
3536// Destructor
3537CbcModel::~CbcModel ()
3538{
3539  if (defaultHandler_) {
3540    delete handler_;
3541    handler_ = NULL;
3542  }
3543  delete tree_;
3544  tree_=NULL;
3545  if (ourSolver_) delete solver_;
3546  gutsOfDestructor();
3547  delete eventHandler_ ;
3548  eventHandler_ = NULL ;
3549}
3550// Clears out as much as possible (except solver)
3551void 
3552CbcModel::gutsOfDestructor()
3553{
3554  delete referenceSolver_;
3555  referenceSolver_=NULL;
3556  int i;
3557  for (i=0;i<numberCutGenerators_;i++) {
3558    delete generator_[i];
3559    delete virginGenerator_[i];
3560  }
3561  delete [] generator_;
3562  delete [] virginGenerator_;
3563  generator_=NULL;
3564  virginGenerator_=NULL;
3565  for (i=0;i<numberHeuristics_;i++)
3566    delete heuristic_[i];
3567  delete [] heuristic_;
3568  heuristic_=NULL;
3569  delete nodeCompare_;
3570  nodeCompare_=NULL;
3571  delete problemFeasibility_;
3572  problemFeasibility_=NULL;
3573  delete [] originalColumns_;
3574  originalColumns_=NULL;
3575  delete strategy_;
3576  gutsOfDestructor2();
3577}
3578// Clears out enough to reset CbcModel
3579void 
3580CbcModel::gutsOfDestructor2()
3581{
3582  delete [] integerInfo_;
3583  integerInfo_=NULL;
3584  delete [] integerVariable_;
3585  integerVariable_=NULL;
3586  int i;
3587  if (ownObjects_) {
3588    for (i=0;i<numberObjects_;i++)
3589      delete object_[i];
3590    delete [] object_;
3591  }
3592  ownObjects_=true;
3593  object_=NULL;
3594  numberIntegers_=0;
3595  numberObjects_=0;
3596  delete emptyWarmStart_ ;
3597  emptyWarmStart_ =NULL;
3598  delete continuousSolver_;
3599  continuousSolver_=NULL;
3600  delete [] bestSolution_;
3601  bestSolution_=NULL;
3602  delete [] currentSolution_;
3603  currentSolution_=NULL;
3604  delete [] continuousSolution_;
3605  continuousSolution_=NULL;
3606  solverCharacteristics_=NULL;
3607  delete [] usedInSolution_;
3608  usedInSolution_ = NULL;
3609  testSolution_=NULL;
3610  lastHeuristic_ = NULL;
3611  delete [] addedCuts_;
3612  addedCuts_=NULL;
3613  nextRowCut_ = NULL;
3614  currentNode_ = NULL;
3615  delete [] walkback_;
3616  walkback_=NULL;
3617  delete [] whichGenerator_;
3618  whichGenerator_ = NULL;
3619  for (i=0;i<maximumStatistics_;i++)
3620    delete statistics_[i];
3621  delete [] statistics_;
3622  statistics_=NULL;
3623  maximumDepthActual_ = 0;
3624  numberDJFixed_ =0.0;
3625  delete probingInfo_;
3626  probingInfo_ = NULL;
3627  maximumStatistics_=0;
3628  delete [] analyzeResults_;
3629  analyzeResults_=NULL;
3630  // Below here is whatever consensus is
3631  ourSolver_=true;
3632  bestObjective_=COIN_DBL_MAX;
3633  bestPossibleObjective_=COIN_DBL_MAX;
3634  sumChangeObjective1_=0.0;
3635  sumChangeObjective2_=0.0;
3636  numberSolutions_=0;
3637  stateOfSearch_=0;
3638  delete [] hotstartSolution_;
3639  hotstartSolution_=NULL;
3640  delete [] hotstartPriorities_;
3641  hotstartPriorities_=NULL;
3642  numberHeuristicSolutions_=0;
3643  numberNodes_=0;
3644  numberNodes2_=0;
3645  numberIterations_=0;
3646  status_=-1;
3647  secondaryStatus_=-1;
3648  maximumNumberCuts_=0;
3649  phase_=0;
3650  currentNumberCuts_=0;
3651  maximumDepth_=0;
3652  nextRowCut_=NULL;
3653  currentNode_=NULL;
3654  delete branchingMethod_;
3655  branchingMethod_=NULL;
3656  delete cutModifier_;
3657  cutModifier_=NULL;
3658  // clear out tree
3659  if (tree_&&tree_->size())
3660    tree_->cleanTree(this, -1.0e100,bestPossibleObjective_) ;
3661  subTreeModel_=NULL;
3662  numberStoppedSubTrees_=0;
3663  numberInfeasibleNodes_=0;
3664  numberGlobalViolations_=0;
3665  continuousObjective_=0.0;
3666  originalContinuousObjective_=0.0;
3667  continuousInfeasibilities_=0;
3668  numberFixedAtRoot_=0;
3669  numberFixedNow_=0;
3670  stoppedOnGap_=false;
3671  eventHappened_=false;
3672  numberLongStrong_=0;
3673  numberOldActiveCuts_=0;
3674  numberNewCuts_=0;
3675  searchStrategy_=-1;
3676  numberStrongIterations_=0;
3677  // Parameters which need to be reset
3678  dblParam_[CbcCutoffIncrement] = 1e-5;
3679  dblParam_[CbcCurrentCutoff] = 1.0e100;
3680  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
3681  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
3682}
3683// Save a copy of the current solver so can be reset to
3684void 
3685CbcModel::saveReferenceSolver()
3686{
3687  delete referenceSolver_;
3688  referenceSolver_= solver_->clone();
3689}
3690
3691// Uses a copy of reference solver to be current solver
3692void 
3693CbcModel::resetToReferenceSolver()
3694{
3695  delete solver_;
3696  solver_ = referenceSolver_->clone();
3697  // clear many things
3698  gutsOfDestructor2();
3699  // Reset cutoff
3700  // Solvers know about direction
3701  double direction = solver_->getObjSense();
3702  double value;
3703  solver_->getDblParam(OsiDualObjectiveLimit,value); 
3704  setCutoff(value*direction);
3705}
3706
3707// Are there a numerical difficulties?
3708bool 
3709CbcModel::isAbandoned() const
3710{
3711  return status_ == 2;
3712}
3713// Is optimality proven?
3714bool 
3715CbcModel::isProvenOptimal() const
3716{
3717  if (!status_ && bestObjective_<1.0e30)
3718    return true;
3719  else
3720    return false;
3721}
3722// Is  infeasiblity proven (or none better than cutoff)?
3723bool 
3724CbcModel::isProvenInfeasible() const
3725{
3726  if (!status_ && bestObjective_>=1.0e30)
3727    return true;
3728  else
3729    return false;
3730}
3731// Was continuous solution unbounded
3732bool 
3733CbcModel::isContinuousUnbounded() const
3734{
3735  if (!status_ && secondaryStatus_==7)
3736    return true;
3737  else
3738    return false;
3739}
3740// Was continuous solution unbounded
3741bool 
3742CbcModel::isProvenDualInfeasible() const
3743{
3744  if (!status_ && secondaryStatus_==7)
3745    return true;
3746  else
3747    return false;
3748}
3749// Node limit reached?
3750bool 
3751CbcModel::isNodeLimitReached() const
3752{
3753  return numberNodes_ >= intParam_[CbcMaxNumNode];
3754}
3755// Time limit reached?
3756bool 
3757CbcModel::isSecondsLimitReached() const
3758{
3759  if (status_==1&&secondaryStatus_==4)
3760    return true;
3761  else
3762    return false;
3763}
3764// Solution limit reached?
3765bool 
3766CbcModel::isSolutionLimitReached() const
3767{
3768  return numberSolutions_ >= intParam_[CbcMaxNumSol];
3769}
3770// Set language
3771void 
3772CbcModel::newLanguage(CoinMessages::Language language)
3773{
3774  messages_ = CbcMessage(language);
3775}
3776void 
3777CbcModel::setNumberStrong(int number)
3778{
3779  if (number<0)
3780    numberStrong_=0;
3781   else
3782    numberStrong_=number;
3783}
3784void 
3785CbcModel::setNumberBeforeTrust(int number)
3786{
3787  if (number<-1) {
3788    numberBeforeTrust_=0;
3789  } else {
3790    numberBeforeTrust_=number;
3791    //numberStrong_ = CoinMax(numberStrong_,1);
3792  }
3793}
3794void 
3795CbcModel::setNumberPenalties(int number)
3796{
3797  if (number<=0) {
3798    numberPenalties_=0;
3799  } else {
3800    numberPenalties_=number;
3801  }
3802}
3803void 
3804CbcModel::setPenaltyScaleFactor(double value)
3805{
3806  if (value<=0) {
3807    penaltyScaleFactor_=3.0;
3808  } else {
3809    penaltyScaleFactor_=value;
3810  }
3811}
3812void 
3813CbcModel::setHowOftenGlobalScan(int number)
3814{
3815  if (number<-1)
3816    howOftenGlobalScan_=0;
3817   else
3818    howOftenGlobalScan_=number;
3819}
3820
3821// Add one generator
3822void 
3823CbcModel::addCutGenerator(CglCutGenerator * generator,
3824                          int howOften, const char * name,
3825                          bool normal, bool atSolution,
3826                          bool whenInfeasible,int howOftenInSub,
3827                          int whatDepth, int whatDepthInSub)
3828{
3829  CbcCutGenerator ** temp = generator_;
3830  generator_ = new CbcCutGenerator * [numberCutGenerators_+1];
3831  memcpy(generator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
3832  delete[] temp ;
3833  generator_[numberCutGenerators_]= 
3834    new CbcCutGenerator(this,generator, howOften, name,
3835                        normal,atSolution,whenInfeasible,howOftenInSub,
3836                        whatDepth, whatDepthInSub);
3837  // and before any cahnges
3838  temp = virginGenerator_;
3839  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_+1];
3840  memcpy(virginGenerator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
3841  delete[] temp ;
3842  virginGenerator_[numberCutGenerators_++]= 
3843    new CbcCutGenerator(this,generator, howOften, name,
3844                        normal,atSolution,whenInfeasible,howOftenInSub,
3845                        whatDepth, whatDepthInSub);
3846                                                         
3847}
3848// Add one heuristic
3849void 
3850CbcModel::addHeuristic(CbcHeuristic * generator, const char *name)
3851{
3852  CbcHeuristic ** temp = heuristic_;
3853  heuristic_ = new CbcHeuristic * [numberHeuristics_+1];
3854  memcpy(heuristic_,temp,numberHeuristics_*sizeof(CbcHeuristic *));
3855  delete [] temp;
3856  heuristic_[numberHeuristics_]=generator->clone();
3857  if (name)
3858  { heuristic_[numberHeuristics_]->setHeuristicName(name) ; }
3859  numberHeuristics_++ ;
3860}
3861
3862/*
3863  The last subproblem handled by the solver is not necessarily related to the
3864  one being recreated, so the first action is to remove all cuts from the
3865  constraint system.  Next, traverse the tree from node to the root to
3866  determine the basis size required for this subproblem and create an empty
3867  basis with the right capacity.  Finally, traverse the tree from root to
3868  node, adjusting bounds in the constraint system, adjusting the basis, and
3869  collecting the cuts that must be added to the constraint system.
3870  applyToModel does the heavy lifting.
3871
3872  addCuts1 is used in contexts where all that's desired is the list of cuts:
3873  the node is already fathomed, and we're collecting cuts so that we can
3874  adjust reference counts as we prune nodes. Arguably the two functions
3875  should be separated. The culprit is applyToModel, which performs cut
3876  collection and model adjustment.
3877
3878  Certainly in the contexts where all we need is a list of cuts, there's no
3879  point in passing in a valid basis --- an empty basis will do just fine.
3880*/
3881void CbcModel::addCuts1 (CbcNode * node, CoinWarmStartBasis *&lastws)
3882{ int i;
3883  int nNode=0;
3884  int numberColumns = getNumCols();
3885  CbcNodeInfo * nodeInfo = node->nodeInfo();
3886
3887/*
3888  Remove all cuts from the constraint system.
3889  (original comment includes ``see note below for later efficiency'', but
3890  the reference isn't clear to me).
3891*/
3892  int currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_;
3893  int *which = new int[currentNumberCuts];
3894  for (i = 0 ; i < currentNumberCuts ; i++)
3895    which[i] = i+numberRowsAtContinuous_;
3896  solver_->deleteRows(currentNumberCuts,which);
3897  delete [] which;
3898/*
3899  Accumulate the path from node to the root in walkback_, and accumulate a
3900  cut count in currentNumberCuts.
3901
3902  original comment: when working then just unwind until where new node joins
3903  old node (for cuts?)
3904*/
3905  currentNumberCuts = 0;
3906  while (nodeInfo) {
3907    //printf("nNode = %d, nodeInfo = %x\n",nNode,nodeInfo);
3908    walkback_[nNode++]=nodeInfo;
3909    currentNumberCuts += nodeInfo->numberCuts() ;
3910    nodeInfo = nodeInfo->parent() ;
3911    if (nNode==maximumDepth_) {
3912      maximumDepth_ *= 2;
3913      CbcNodeInfo ** temp = new CbcNodeInfo * [maximumDepth_];
3914      for (i=0;i<nNode;i++) 
3915        temp[i] = walkback_[i];
3916      delete [] walkback_;
3917      walkback_ = temp;
3918    }
3919  }
3920/*
3921  Create an empty basis with sufficient capacity for the constraint system
3922  we'll construct: original system plus cuts. Make sure we have capacity to
3923  record those cuts in addedCuts_.
3924
3925  The method of adjusting the basis at a FullNodeInfo object (the root, for
3926  example) is to use a copy constructor to duplicate the basis held in the
3927  nodeInfo, then resize it and return the new basis object. Guaranteed,
3928  lastws will point to a different basis when it returns. We pass in a basis
3929  because we need the parameter to return the allocated basis, and it's an
3930  easy way to pass in the size. But we take a hit for memory allocation.
3931*/
3932  currentNumberCuts_=currentNumberCuts;
3933  if (currentNumberCuts > maximumNumberCuts_) {
3934    maximumNumberCuts_ = currentNumberCuts;
3935    delete [] addedCuts_;
3936    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3937  }
3938  lastws->setSize(numberColumns,numberRowsAtContinuous_+currentNumberCuts);
3939/*
3940  This last bit of code traverses the path collected in walkback_ from the
3941  root back to node. At the end of the loop,
3942   * lastws will be an appropriate basis for node;
3943   * variable bounds in the constraint system will be set to be correct for
3944     node; and
3945   * addedCuts_ will be set to a list of cuts that need to be added to the
3946     constraint system at node.
3947  applyToModel does all the heavy lifting.
3948*/
3949  currentNumberCuts=0;
3950  while (nNode) {
3951    --nNode;
3952    walkback_[nNode]->applyToModel(this,lastws,addedCuts_,currentNumberCuts);
3953  }
3954}
3955
3956/*
3957  adjustCuts might be a better name: If the node is feasible, we sift through
3958  the cuts collected by addCuts1, add the ones that are tight and omit the
3959  ones that are loose. If the node is infeasible, we just adjust the
3960  reference counts to reflect that we're about to prune this node and its
3961  descendants.
3962*/
3963int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws,bool canFix)
3964{
3965/*
3966  addCuts1 performs step 1 of restoring the subproblem at this node; see the
3967  comments there.
3968*/
3969  addCuts1(node,lastws);
3970  int i;
3971  int numberColumns = getNumCols();
3972  CbcNodeInfo * nodeInfo = node->nodeInfo();
3973  double cutoff = getCutoff() ;
3974  int currentNumberCuts=currentNumberCuts_;
3975  if (canFix) {
3976    bool feasible=true;
3977    const double *lower = solver_->getColLower() ;
3978    const double *upper = solver_->getColUpper() ;
3979    double * newLower = analyzeResults_;
3980    double * objLower = newLower+numberIntegers_;
3981    double * newUpper = objLower+numberIntegers_;
3982    double * objUpper = newUpper+numberIntegers_;
3983    int n=0;
3984    for (i=0;i<numberIntegers_;i++) {
3985      int iColumn = integerVariable_[i];
3986      bool changed=false;
3987      double lo = 0.0;
3988      double up = 0.0;
3989      if (objLower[i]>cutoff) {
3990        lo = lower[iColumn];
3991        up = upper[iColumn];
3992        if (lo<newLower[i]) {
3993          lo = newLower[i];
3994          solver_->setColLower(iColumn,lo);
3995          changed=true;
3996          n++;
3997        }
3998        if (objUpper[i]>cutoff) {
3999          if (up>newUpper[i]) {
4000            up = newUpper[i];
4001            solver_->setColUpper(iColumn,up);
4002            changed=true;
4003            n++;
4004          }
4005        }
4006      } else if (objUpper[i]>cutoff) {
4007        lo = lower[iColumn];
4008        up = upper[iColumn];
4009        if (up>newUpper[i]) {
4010          up = newUpper[i];
4011          solver_->setColUpper(iColumn,up);
4012          changed=true;
4013          n++;
4014        }
4015      }
4016      if (changed&&lo>up) {
4017        feasible=false;
4018        break;
4019      }
4020    }
4021    if (!feasible) {
4022      printf("analysis says node infeas\n");
4023      cutoff=-COIN_DBL_MAX;
4024    }
4025  }
4026/*
4027  If the node can't be fathomed by bound, reinstall tight cuts in the
4028  constraint system. Even if there are no cuts, we'll want to set the
4029  reconstructed basis in the solver.
4030*/
4031  if (node->objectiveValue() < cutoff||numberThreads_)
4032  { 
4033#   ifdef CBC_CHECK_BASIS
4034    printf("addCuts: expanded basis; rows %d+%d\n",
4035           numberRowsAtContinuous_,currentNumberCuts);
4036    lastws->print();
4037#   endif
4038/*
4039  Adjust the basis and constraint system so that we retain only active cuts.
4040  There are three steps:
4041    1) Scan the basis. Sort the cuts into effective cuts to be kept and
4042       loose cuts to be dropped.
4043    2) Drop the loose cuts and resize the basis to fit.
4044    3) Install the tight cuts in the constraint system (applyRowCuts) and
4045       and install the basis (setWarmStart).
4046  Use of compressRows conveys we're compressing the basis and not just
4047  tweaking the artificialStatus_ array.
4048*/
4049    if (currentNumberCuts > 0) {
4050      int numberToAdd = 0;
4051      const OsiRowCut **addCuts;
4052      int numberToDrop = 0 ;
4053      int *cutsToDrop ;
4054      addCuts = new const OsiRowCut* [currentNumberCuts];
4055      cutsToDrop = new int[currentNumberCuts] ;
4056      assert (currentNumberCuts+numberRowsAtContinuous_<=lastws->getNumArtificial());
4057      for (i=0;i<currentNumberCuts;i++) {
4058        CoinWarmStartBasis::Status status = 
4059          lastws->getArtifStatus(i+numberRowsAtContinuous_);
4060        if (addedCuts_[i] &&
4061            (status != CoinWarmStartBasis::basic ||
4062             addedCuts_[i]->effectiveness()==COIN_DBL_MAX)) {
4063#         ifdef CHECK_CUT_COUNTS
4064          printf("Using cut %d %x as row %d\n",i,addedCuts_[i],
4065                 numberRowsAtContinuous_+numberToAdd);
4066#         endif
4067          addCuts[numberToAdd++] = new OsiRowCut(*addedCuts_[i]);
4068        } else {
4069#         ifdef CHECK_CUT_COUNTS
4070          printf("Dropping cut %d %x\n",i,addedCuts_[i]);
4071#         endif
4072          addedCuts_[i]=NULL;
4073          cutsToDrop[numberToDrop++] = numberRowsAtContinuous_+i ;
4074        }
4075      }
4076      int numberRowsNow=numberRowsAtContinuous_+numberToAdd;
4077      lastws->compressRows(numberToDrop,cutsToDrop) ;
4078      lastws->resize(numberRowsNow,numberColumns);
4079      solver_->applyRowCuts(numberToAdd,addCuts);
4080#     ifdef CBC_CHECK_BASIS
4081      printf("addCuts: stripped basis; rows %d + %d\n",
4082             numberRowsAtContinuous_,numberToAdd);
4083      lastws->print();
4084#     endif
4085      for (i=0;i<numberToAdd;i++)
4086        delete addCuts[i];
4087      delete [] addCuts;
4088      delete [] cutsToDrop ;
4089    }
4090/*
4091  Set the basis in the solver.
4092*/
4093    solver_->setWarmStart(lastws);
4094#if 0
4095    if ((numberNodes_%printFrequency_)==0) {
4096      printf("Objective %g, depth %d, unsatisfied %d\n",
4097             node->objectiveValue(),
4098             node->depth(),node->numberUnsatisfied());
4099    }
4100#endif
4101/*
4102  Clean up and we're out of here.
4103*/
4104    numberNodes_++;
4105    return 0;
4106  } 
4107/*
4108  This node has been fathomed by bound as we try to revive it out of the live
4109  set. Adjust the cut reference counts to reflect that we no longer need to
4110  explore the remaining branch arms, hence they will no longer reference any
4111  cuts. Cuts whose reference count falls to zero are deleted. 
4112*/
4113  else
4114  { int i;
4115    if (currentNumberCuts) {
4116      lockThread();
4117      int numberLeft = nodeInfo->numberBranchesLeft();
4118      for (i = 0 ; i < currentNumberCuts ; i++)
4119        { if (addedCuts_[i])
4120          { if (!addedCuts_[i]->decrement(numberLeft))
4121            { delete addedCuts_[i];
4122            addedCuts_[i] = NULL; } } }
4123      unlockThread();
4124    }
4125    return 1 ; }
4126}
4127
4128
4129/*
4130  Perform reduced cost fixing on integer variables.
4131
4132  The variables in question are already nonbasic at bound. We're just nailing
4133  down the current situation.
4134*/
4135
4136int CbcModel::reducedCostFix ()
4137
4138{
4139  if(!solverCharacteristics_->reducedCostsAccurate())
4140    return 0; //NLP
4141  double cutoff = getCutoff() ;
4142  double direction = solver_->getObjSense() ;
4143  double gap = cutoff - solver_->getObjValue()*direction ;
4144  double tolerance;
4145  solver_->getDblParam(OsiDualTolerance,tolerance) ;
4146  if (gap<=0.0)
4147    return 0;
4148  gap += 100.0*tolerance;
4149  double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4150
4151  const double *lower = solver_->getColLower() ;
4152  const double *upper = solver_->getColUpper() ;
4153  const double *solution = solver_->getColSolution() ;
4154  const double *reducedCost = solver_->getReducedCost() ;
4155
4156  int numberFixed = 0 ;
4157
4158# ifdef COIN_HAS_CLP
4159  OsiClpSolverInterface * clpSolver
4160    = dynamic_cast<OsiClpSolverInterface *> (solver_);
4161  ClpSimplex * clpSimplex=NULL;
4162  if (clpSolver) 
4163    clpSimplex = clpSolver->getModelPtr();
4164# endif
4165  for (int i = 0 ; i < numberIntegers_ ; i++)
4166  { int iColumn = integerVariable_[i] ;
4167    double djValue = direction*reducedCost[iColumn] ;
4168    if (upper[iColumn]-lower[iColumn] > integerTolerance)
4169    { if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap)
4170      { solver_->setColUpper(iColumn,lower[iColumn]) ;
4171#ifdef COIN_HAS_CLP
4172      // may just have been fixed before
4173      if (clpSimplex)
4174        assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
4175               clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4176#endif
4177      numberFixed++ ; }
4178      else
4179      if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap)
4180      { solver_->setColLower(iColumn,upper[iColumn]) ;
4181#ifdef COIN_HAS_CLP
4182      // may just have been fixed before
4183      if (clpSimplex)
4184        assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
4185               clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4186#endif
4187      numberFixed++ ; } } }
4188  numberDJFixed_ += numberFixed;
4189 
4190  return numberFixed; }
4191
4192// Collect coding to replace whichGenerator
4193void
4194CbcModel::resizeWhichGenerator(int numberNow, int numberAfter)
4195{
4196  if (numberAfter > maximumWhich_) {
4197    maximumWhich_ = CoinMax(maximumWhich_*2+100,numberAfter) ;
4198    int * temp = new int[2*maximumWhich_] ;
4199    memcpy(temp,whichGenerator_,numberNow*sizeof(int)) ;
4200    delete [] whichGenerator_ ;
4201    whichGenerator_ = temp ;
4202    memset(whichGenerator_+numberNow,0,(maximumWhich_-numberNow)*sizeof(int));
4203  }
4204}
4205
4206/** Solve the model using cuts
4207
4208  This version takes off redundant cuts from node.
4209  Returns true if feasible.
4210
4211  \todo
4212  Why do I need to resolve the problem? What has been done between the last
4213  relaxation and calling solveWithCuts?
4214
4215  If numberTries == 0 then user did not want any cuts.
4216*/
4217
4218bool 
4219CbcModel::solveWithCuts (OsiCuts &cuts, int numberTries, CbcNode *node)
4220/*
4221  Parameters:
4222    numberTries: (i) the maximum number of iterations for this round of cut
4223                     generation; if negative then we don't mind if drop is tiny.
4224   
4225    cuts:       (o) all cuts generated in this round of cut generation
4226
4227    node: (i)     So we can update dynamic pseudo costs
4228*/
4229                       
4230
4231{
4232# ifdef COIN_HAS_CLP
4233  OsiClpSolverInterface * clpSolver
4234    = dynamic_cast<OsiClpSolverInterface *> (solver_);
4235  int saveClpOptions=0;
4236  if (clpSolver) 
4237    saveClpOptions = clpSolver->specialOptions();
4238# endif
4239#ifdef CBC_THREAD
4240  CbcModel ** threadModel = NULL;
4241  pthread_t * threadId = NULL;
4242  pthread_cond_t condition_main;
4243  pthread_mutex_t condition_mutex;
4244  pthread_mutex_t * mutex2 = NULL;
4245  pthread_cond_t * condition2 = NULL;
4246  threadStruct * threadInfo = NULL;
4247  void * saveMutex = NULL;
4248  if (numberThreads_&&(threadMode_&2)!=0&&!numberNodes_) {
4249    threadId = new pthread_t [numberThreads_];
4250    pthread_cond_init(&condition_main,NULL);
4251    pthread_mutex_init(&condition_mutex,NULL);
4252    threadModel = new CbcModel * [numberThreads_];
4253    threadInfo = new threadStruct [numberThreads_+1];
4254    mutex2 = new pthread_mutex_t [numberThreads_];
4255    condition2 = new pthread_cond_t [numberThreads_];
4256    saveMutex = mutex_;
4257    for (int i=0;i<numberThreads_;i++) {
4258      pthread_mutex_init(mutex2+i,NULL);
4259      pthread_cond_init(condition2+i,NULL);
4260      threadId[i]=0;
4261      threadModel[i]=new CbcModel;
4262      threadModel[i]->generator_ = new CbcCutGenerator * [1];
4263      delete threadModel[i]->solver_;
4264      threadModel[i]->solver_=NULL;
4265      threadModel[i]->numberThreads_=numberThreads_;
4266      mutex_ = (void *) (threadInfo+i);
4267      threadInfo[i].thisModel=(CbcModel *) threadModel[i];
4268      threadInfo[i].baseModel=this;
4269      threadInfo[i].threadIdOfBase=pthread_self();
4270      threadInfo[i].mutex2=mutex2+i;
4271      threadInfo[i].condition2=condition2+i;
4272      threadInfo[i].returnCode=-1;
4273      pthread_create(threadId+i,NULL,doCutsThread,threadInfo+i);
4274    }
4275    // Do a partial one for base model
4276    threadInfo[numberThreads_].baseModel=this;
4277    mutex_ = (void *) (threadInfo+numberThreads_);
4278    threadInfo[numberThreads_].condition2=&condition_main;
4279    threadInfo[numberThreads_].mutex2=&condition_mutex;
4280  }
4281#endif
4282  bool feasible = true ;
4283  int lastNumberCuts = 0 ;
4284  double lastObjective = -1.0e100 ;
4285  int violated = 0 ;
4286  int numberRowsAtStart = solver_->getNumRows() ;
4287  //printf("solver had %d rows\n",numberRowsAtStart);
4288  int numberColumns = solver_->getNumCols() ;
4289  CoinBigIndex numberElementsAtStart = solver_->getNumElements();
4290
4291  numberOldActiveCuts_ = numberRowsAtStart-numberRowsAtContinuous_ ;
4292  numberNewCuts_ = 0 ;
4293
4294  bool onOptimalPath = false ;
4295  const OsiRowCutDebugger *debugger = NULL;
4296  if ((specialOptions_&1)!=0) {
4297    /*
4298      See OsiRowCutDebugger for details. In a nutshell, make sure that current
4299      variable values do not conflict with a known optimal solution. (Obviously
4300      this can be fooled when there are multiple solutions.)
4301    */
4302    debugger = solver_->getRowCutDebugger() ;
4303    if (debugger) 
4304      onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
4305  }
4306  OsiCuts slackCuts;
4307/*
4308  Resolve the problem. If we've lost feasibility, might as well bail out right
4309  after the debug stuff. The resolve will also refresh cached copies of the
4310  solver solution (cbcColLower_, ...) held by CbcModel.
4311*/
4312  double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
4313  if (node)
4314    objectiveValue= node->objectiveValue();
4315  int returnCode = resolve(node ? node->nodeInfo() : NULL,1);
4316#ifdef COIN_DEVELOP
4317  //if (!solver_->getIterationCount()&&solver_->isProvenOptimal())
4318  //printf("zero iterations on first solve of branch\n");
4319#endif
4320  if (node&&node->nodeInfo()&&!node->nodeInfo()->numberBranchesLeft())
4321    node->nodeInfo()->allBranchesGone(); // can clean up
4322  feasible = returnCode  != 0 ;
4323  if (returnCode<0)
4324    numberTries=0;
4325  if (problemFeasibility_->feasible(this,0)<0) {
4326    feasible=false; // pretend infeasible
4327  }
4328 
4329  // Update branching information if wanted
4330  if(node &&branchingMethod_)
4331    branchingMethod_->updateInformation(solver_,node);
4332
4333#ifdef CBC_DEBUG
4334  if (feasible)
4335  { printf("Obj value %g (%s) %d rows\n",solver_->getObjValue(),
4336           (solver_->isProvenOptimal())?"proven":"unproven",
4337           solver_->getNumRows()) ; }
4338 
4339  else
4340  { printf("Infeasible %d rows\n",solver_->getNumRows()) ; }
4341#endif
4342  if ((specialOptions_&1)!=0) {
4343/*
4344  If the RowCutDebugger said we were compatible with the optimal solution,
4345  and now we're suddenly infeasible, we might be confused. Then again, we
4346  may have fathomed by bound, heading for a rediscovery of an optimal solution.
4347*/
4348    if (onOptimalPath && !solver_->isDualObjectiveLimitReached()) {
4349      if (!feasible) {
4350        solver_->writeMps("infeas");
4351        CoinWarmStartBasis *slack =
4352          dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
4353        solver_->setWarmStart(slack);
4354        delete slack ;
4355        solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
4356        solver_->initialSolve();
4357      }
4358      assert(feasible) ;
4359    }
4360  }
4361
4362  if (!feasible) {
4363    numberInfeasibleNodes_++;
4364# ifdef COIN_HAS_CLP
4365    if (clpSolver) 
4366    clpSolver->setSpecialOptions(saveClpOptions);
4367# endif
4368    return (false) ;
4369  }
4370  sumChangeObjective1_ += solver_->getObjValue()*solver_->getObjSense()
4371    - objectiveValue ;
4372  if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
4373    numberTries=0; // exit
4374  //if ((numberNodes_%100)==0)
4375  //printf("XXa sum obj changed by %g\n",sumChangeObjective1_);
4376  objectiveValue = solver_->getObjValue()*solver_->getObjSense();
4377  // Return at once if numberTries zero
4378  if (!numberTries) {
4379    cuts=OsiCuts();
4380    numberNewCuts_=0;
4381# ifdef COIN_HAS_CLP
4382    if (clpSolver) 
4383    clpSolver->setSpecialOptions(saveClpOptions);
4384# endif
4385    return true;
4386  }
4387/*
4388  Do reduced cost fixing.
4389*/
4390  reducedCostFix() ;
4391/*
4392  Set up for at most numberTries rounds of cut generation. If numberTries is
4393  negative, we'll ignore the minimumDrop_ cutoff and keep generating cuts for
4394  the specified number of rounds.
4395*/
4396  double minimumDrop = minimumDrop_ ;
4397  if (numberTries<0)
4398  { numberTries = -numberTries ;
4399    minimumDrop = -1.0 ; }
4400/*
4401  Is it time to scan the cuts in order to remove redundant cuts? If so, set
4402  up to do it.
4403*/
4404# define SCANCUTS 100 
4405  int *countColumnCuts = NULL ;
4406  // Always accumulate row cut counts
4407  int * countRowCuts =new int[numberCutGenerators_+numberHeuristics_] ;
4408  memset(countRowCuts,0,
4409         (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ;
4410  bool fullScan = false ;
4411  if ((numberNodes_%SCANCUTS) == 0)
4412  { fullScan = true ;
4413    countColumnCuts = new int[numberCutGenerators_+numberHeuristics_] ;
4414    memset(countColumnCuts,0,
4415           (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ; }
4416
4417  double direction = solver_->getObjSense() ;
4418  double startObjective = solver_->getObjValue()*direction ;
4419
4420  currentPassNumber_ = 0 ;
4421  double primalTolerance = 1.0e-7 ;
4422  // We may need to keep going on
4423  bool keepGoing=false;
4424/*
4425  Begin cut generation loop. Cuts generated during each iteration are
4426  collected in theseCuts. The loop can be divided into four phases:
4427   1) Prep: Fix variables using reduced cost. In the first iteration only,
4428      consider scanning globalCuts_ and activating any applicable cuts.
4429   2) Cut Generation: Call each generator and heuristic registered in the
4430      generator_ and heuristic_ arrays. Newly generated global cuts are
4431      copied to globalCuts_ at this time.
4432   3) Cut Installation and Reoptimisation: Install column and row cuts in
4433      the solver. Copy row cuts to cuts (parameter). Reoptimise.
4434   4) Cut Purging: takeOffCuts() removes inactive cuts from the solver, and
4435      does the necessary bookkeeping in the model.
4436*/
4437  do
4438  { currentPassNumber_++ ;
4439    numberTries-- ;
4440    if (numberTries<0&&keepGoing) {
4441      // switch off all normal ones
4442      for (int i = 0;i<numberCutGenerators_;i++) {
4443        if (!generator_[i]->mustCallAgain())
4444          generator_[i]->setSwitchedOff(true);
4445      }
4446    }
4447    keepGoing=false;
4448    OsiCuts theseCuts ;
4449/*
4450  Scan previously generated global column and row cuts to see if any are
4451  useful.
4452*/
4453    int numberViolated=0;
4454    if (currentPassNumber_ == 1 && howOftenGlobalScan_ > 0 &&
4455        (numberNodes_%howOftenGlobalScan_) == 0)
4456    { int numberCuts = globalCuts_.sizeColCuts() ;
4457      int i;
4458      // possibly extend whichGenerator
4459      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
4460      for ( i = 0 ; i < numberCuts ; i++)
4461      { OsiColCut *thisCut = globalCuts_.colCutPtr(i) ;
4462        if (thisCut->violated(cbcColSolution_)>primalTolerance) {
4463          printf("Global cut added - violation %g\n",
4464                 thisCut->violated(cbcColSolution_)) ;
4465          whichGenerator_[numberViolated++]=-1;
4466#ifndef GLOBAL_CUTS_JUST_POINTERS
4467          theseCuts.insert(*thisCut) ;
4468#else
4469          theseCuts.insert(thisCut) ;
4470#endif
4471        }
4472      }
4473      numberCuts = globalCuts_.sizeRowCuts() ;
4474      // possibly extend whichGenerator
4475      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
4476      for ( i = 0;i<numberCuts;i++) {
4477        OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ;
4478        if (thisCut->violated(cbcColSolution_)>primalTolerance) {
4479          //printf("Global cut added - violation %g\n",
4480          // thisCut->violated(cbcColSolution_)) ;
4481          whichGenerator_[numberViolated++]=-1;
4482#ifndef GLOBAL_CUTS_JUST_POINTERS
4483          theseCuts.insert(*thisCut) ;
4484#else
4485          theseCuts.insert(thisCut) ;
4486#endif
4487        }
4488      }
4489      numberGlobalViolations_+=numberViolated;
4490    }
4491/*
4492  Generate new cuts (global and/or local) and/or apply heuristics.  If
4493  CglProbing is used, then it should be first as it can fix continuous
4494  variables.
4495
4496  At present, CglProbing is the only case where generateCuts will return
4497  true. generateCuts actually modifies variable bounds in the solver when
4498  CglProbing indicates that it can fix a variable. Reoptimisation is required
4499  to take full advantage.
4500
4501  The need to resolve here should only happen after a heuristic solution.
4502  (Note default OSI implementation of optimalBasisIsAvailable always returns
4503  false.)
4504*/
4505    if (solverCharacteristics_->warmStart()&&
4506        !solver_->optimalBasisIsAvailable()) {
4507      //printf("XXXXYY no opt basis\n");
4508      resolve(node ? node->nodeInfo() : NULL,3);
4509    }
4510    if (nextRowCut_) {
4511      // branch was a cut - add it
4512      theseCuts.insert(*nextRowCut_);
4513      if (handler_->logLevel()>1)
4514        nextRowCut_->print();
4515      const OsiRowCut * cut=nextRowCut_;
4516      double lb = cut->lb();
4517      double ub = cut->ub();
4518      int n=cut->row().getNumElements();
4519      const int * column = cut->row().getIndices();
4520      const double * element = cut->row().getElements();
4521      double sum=0.0;
4522      for (int i=0;i<n;i++) {
4523        int iColumn = column[i];
4524        double value = element[i];
4525        //if (cbcColSolution_[iColumn]>1.0e-7)
4526        //printf("value of %d is %g\n",iColumn,cbcColSolution_[iColumn]);
4527        sum += value * cbcColSolution_[iColumn];
4528      }
4529      delete nextRowCut_;
4530      nextRowCut_=NULL;
4531      if (handler_->logLevel()>1)
4532        printf("applying branch cut, sum is %g, bounds %g %g\n",sum,lb,ub);
4533      // possibly extend whichGenerator
4534      resizeWhichGenerator(numberViolated, numberViolated+1);
4535      // set whichgenerator (also serves as marker to say don't delete0
4536      whichGenerator_[numberViolated++]=-2;
4537    }
4538
4539    // reset probing info
4540    if (probingInfo_)
4541      probingInfo_->initializeFixing();
4542    int i;
4543    if ((threadMode_&2)==0||numberNodes_) {
4544      for (i = 0;i<numberCutGenerators_;i++) {
4545        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
4546        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
4547        bool generate = generator_[i]->normal();
4548        // skip if not optimal and should be (maybe a cut generator has fixed variables)
4549        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
4550          generate=false;
4551        if (generator_[i]->switchedOff())
4552          generate=false;;
4553        if (generate) {
4554          bool mustResolve = 
4555            generator_[i]->generateCuts(theseCuts,fullScan,solver_,node) ;
4556          int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
4557          if(numberRowCutsBefore < numberRowCutsAfter &&
4558             generator_[i]->mustCallAgain())
4559            keepGoing=true; // say must go round
4560          // Check last cut to see if infeasible
4561          if(numberRowCutsBefore < numberRowCutsAfter) {
4562            const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
4563            if (thisCut->lb()>thisCut->ub()) {
4564              feasible = false; // sub-problem is infeasible
4565              break;
4566            }
4567          }
4568#ifdef CBC_DEBUG
4569          {
4570            int k ;
4571            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
4572              OsiRowCut thisCut = theseCuts.rowCut(k) ;
4573              /* check size of elements.
4574                 We can allow smaller but this helps debug generators as it
4575                 is unsafe to have small elements */
4576              int n=thisCut.row().getNumElements();
4577              const int * column = thisCut.row().getIndices();
4578              const double * element = thisCut.row().getElements();
4579              //assert (n);
4580              for (int i=0;i<n;i++) {
4581                double value = element[i];
4582                assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
4583              }
4584            }
4585          }
4586#endif
4587          if (mustResolve) {
4588            int returncode = resolve(node ? node->nodeInfo() : NULL,2);
4589            feasible = returnCode  != 0 ;
4590            if (returncode<0)
4591              numberTries=0;
4592            if ((specialOptions_&1)!=0) {
4593              debugger = solver_->getRowCutDebugger() ;
4594              if (debugger) 
4595                onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
4596              else
4597                onOptimalPath=false;
4598              if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
4599                assert(feasible) ;
4600            }
4601            if (!feasible)
4602              break ;
4603          }
4604        }
4605        int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
4606        int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
4607       
4608        if ((specialOptions_&1)!=0) {
4609          if (onOptimalPath) {
4610            int k ;
4611            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
4612              OsiRowCut thisCut = theseCuts.rowCut(k) ;
4613              if(debugger->invalidCut(thisCut)) {
4614                solver_->writeMps("badCut");
4615#ifdef NDEBUG
4616                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
4617                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
4618                abort();
4619#endif
4620              }
4621              assert(!debugger->invalidCut(thisCut)) ;
4622            }
4623          }
4624        }
4625/*
4626  The cut generator has done its thing, and maybe it generated some
4627  cuts.  Do a bit of bookkeeping: load
4628  whichGenerator[i] with the index of the generator responsible for a cut,
4629  and place cuts flagged as global in the global cut pool for the model.
4630
4631  lastNumberCuts is the sum of cuts added in previous iterations; it's the
4632  offset to the proper starting position in whichGenerator.
4633*/
4634        int numberBefore =
4635          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
4636        int numberAfter =
4637          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
4638        // possibly extend whichGenerator
4639        resizeWhichGenerator(numberBefore, numberAfter);
4640        int j ;
4641        if (fullScan) {
4642          // counts
4643          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
4644        }
4645        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
4646       
4647        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
4648          whichGenerator_[numberBefore++] = i ;
4649          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
4650          if (thisCut->lb()>thisCut->ub())
4651            violated=-2; // sub-problem is infeasible
4652          if (thisCut->globallyValid()) {
4653            // add to global list
4654            OsiRowCut newCut(*thisCut);
4655            newCut.setGloballyValid(2);
4656            newCut.mutableRow().setTestForDuplicateIndex(false);
4657            globalCuts_.insert(newCut) ;
4658          }
4659        }
4660        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
4661          whichGenerator_[numberBefore++] = i ;
4662          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
4663          if (thisCut->globallyValid()) {
4664            // add to global list
4665            OsiColCut newCut(*thisCut);
4666            newCut.setGloballyValid(2);
4667            globalCuts_.insert(newCut) ;
4668          }
4669        }
4670      }
4671      // Add in any violated saved cuts
4672      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
4673        int numberOld = theseCuts.sizeRowCuts();
4674        int numberCuts = slackCuts.sizeRowCuts() ;
4675        int i;
4676        // possibly extend whichGenerator
4677        resizeWhichGenerator(numberOld, numberOld+numberCuts);
4678        for ( i = 0;i<numberCuts;i++) {
4679          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
4680          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
4681            if (messageHandler()->logLevel()>2)
4682              printf("Old cut added - violation %g\n",
4683                     thisCut->violated(cbcColSolution_)) ;
4684            whichGenerator_[numberOld++]=-1;
4685            theseCuts.insert(*thisCut) ;
4686          }
4687        }
4688      }
4689    } else {
4690      // do cuts independently
4691      OsiCuts * eachCuts = new OsiCuts [numberCutGenerators_];;
4692#ifdef CBC_THREAD
4693      if (!threadModel) {
4694#endif
4695        // generate cuts
4696        for (i = 0;i<numberCutGenerators_;i++) {
4697          bool generate = generator_[i]->normal();
4698          // skip if not optimal and should be (maybe a cut generator has fixed variables)
4699          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
4700            generate=false;
4701          if (generator_[i]->switchedOff())
4702            generate=false;;
4703          if (generate) 
4704            generator_[i]->generateCuts(eachCuts[i],fullScan,solver_,node) ;
4705        }
4706#ifdef CBC_THREAD
4707      } else {
4708        for (i=0;i<numberThreads_;i++) {
4709          // set solver here after cloning
4710          threadModel[i]->solver_=solver_->clone();
4711          threadModel[i]->numberNodes_ = (fullScan) ? 1 : 0;
4712        }
4713        // generate cuts
4714        for (i = 0;i<numberCutGenerators_;i++) {
4715          bool generate = generator_[i]->normal();
4716          // skip if not optimal and should be (maybe a cut generator has fixed variables)
4717          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
4718            generate=false;
4719          if (generator_[i]->switchedOff())
4720            generate=false;;
4721          if (generate) { 
4722            bool finished=false;
4723            int iThread=-1;
4724            // see if any available
4725            for (iThread=0;iThread<numberThreads_;iThread++) {
4726              if (threadInfo[iThread].returnCode) {
4727                finished=true;
4728                break;
4729              } else if (threadInfo[iThread].returnCode==0) {
4730#ifdef CONDITION_WAIT
4731                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4732#endif
4733              }
4734            }
4735            while (!finished) {
4736#ifdef CONDITION_WAIT
4737              pthread_mutex_lock(&condition_mutex);
4738              struct timespec absTime;
4739              clock_gettime(CLOCK_REALTIME,&absTime);
4740              absTime.tv_nsec += 1000000; // millisecond
4741              if (absTime.tv_nsec>=1000000000) {
4742                absTime.tv_nsec -= 1000000000;
4743                absTime.tv_sec++;
4744              }
4745              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
4746              pthread_mutex_unlock(&condition_mutex);
4747#else
4748              struct timespec moreTime;
4749              moreTime.tv_nsec = 100;
4750              moreTime.tv_sec = 0 ;
4751              nanosleep(&moreTime,NULL);
4752#endif
4753              for (iThread=0;iThread<numberThreads_;iThread++) {
4754                if (threadInfo[iThread].returnCode>0) {
4755                  finished=true;
4756                  break;
4757                } else if (threadInfo[iThread].returnCode==0) {
4758#ifdef CONDITION_WAIT
4759                  pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4760#endif
4761                }
4762              }
4763            }
4764            assert (iThread<numberThreads_);
4765            assert (threadInfo[iThread].returnCode);
4766            threadModel[iThread]->generator_[0]=generator_[i];
4767            threadModel[iThread]->object_ = (OsiObject **) (eachCuts+i);
4768            // allow to start
4769            threadInfo[iThread].returnCode=0;
4770#ifdef CONDITION_WAIT
4771            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4772#endif
4773          }
4774        }
4775        // wait
4776        for (int iThread=0;iThread<numberThreads_;iThread++) {
4777          if (threadInfo[iThread].returnCode==0) {
4778            bool finished=false;
4779#ifdef CONDITION_WAIT
4780            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4781#endif
4782            while (!finished) {
4783#ifdef CONDITION_WAIT
4784              pthread_mutex_lock(&condition_mutex);
4785              struct timespec absTime;
4786              clock_gettime(CLOCK_REALTIME,&absTime);
4787              absTime.tv_nsec += 1000000; // millisecond
4788              if (absTime.tv_nsec>=1000000000) {
4789                absTime.tv_nsec -= 1000000000;
4790                absTime.tv_sec++;
4791              }
4792              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
4793              pthread_mutex_unlock(&condition_mutex);
4794#else
4795              struct timespec moreTime;
4796              moreTime.tv_nsec = 100;
4797              moreTime.tv_sec = 0;
4798              nanosleep(&moreTime,NULL);
4799#endif
4800              if (threadInfo[iThread].returnCode>0) {
4801                finished=true;
4802                break;
4803              } else if (threadInfo[iThread].returnCode==0) {
4804#ifdef CONDITION_WAIT
4805                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4806#endif
4807              }
4808            }
4809          }
4810          assert (threadInfo[iThread].returnCode);
4811          // say available
4812          threadInfo[iThread].returnCode=-1;
4813          delete threadModel[iThread]->solver_;
4814          threadModel[iThread]->solver_=NULL;
4815        }
4816      }
4817#endif
4818      // Now put together
4819      for (i = 0;i<numberCutGenerators_;i++) {
4820        // add column cuts
4821        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
4822        int numberColumnCuts = eachCuts[i].sizeColCuts();
4823        int numberColumnCutsAfter = numberColumnCutsBefore
4824          + numberColumnCuts;
4825        int j;
4826        for (j=0;j<numberColumnCuts;j++) {
4827          theseCuts.insert(eachCuts[i].colCut(j));
4828        }
4829        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
4830        int numberRowCuts = eachCuts[i].sizeRowCuts();
4831        int numberRowCutsAfter = numberRowCutsBefore
4832          + numberRowCuts;
4833        if (numberRowCuts) {
4834          for (j=0;j<numberRowCuts;j++) {
4835            theseCuts.insert(eachCuts[i].rowCut(j));
4836          }
4837          if (generator_[i]->mustCallAgain())
4838            keepGoing=true; // say must go round
4839          // Check last cut to see if infeasible
4840          const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
4841          if (thisCut->lb()>thisCut->ub()) {
4842            feasible = false; // sub-problem is infeasible
4843            break;
4844          }
4845        }
4846#ifdef CBC_DEBUG
4847        {
4848          int k ;
4849          for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
4850            OsiRowCut thisCut = theseCuts.rowCut(k) ;
4851            /* check size of elements.
4852               We can allow smaller but this helps debug generators as it
4853               is unsafe to have small elements */
4854            int n=thisCut.row().getNumElements();
4855            const int * column = thisCut.row().getIndices();
4856            const double * element = thisCut.row().getElements();
4857            //assert (n);
4858            for (int i=0;i<n;i++) {
4859              double value = element[i];
4860              assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
4861            }
4862          }
4863        }
4864#endif
4865        if ((specialOptions_&1)!=0) {
4866          if (onOptimalPath) {
4867            int k ;
4868            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
4869              OsiRowCut thisCut = theseCuts.rowCut(k) ;
4870              if(debugger->invalidCut(thisCut)) {
4871                solver_->writeMps("badCut");
4872#ifdef NDEBUG
4873                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
4874                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
4875                abort();
4876#endif
4877              }
4878              assert(!debugger->invalidCut(thisCut)) ;
4879            }
4880          }
4881        }
4882/*
4883  The cut generator has done its thing, and maybe it generated some
4884  cuts.  Do a bit of bookkeeping: load
4885  whichGenerator[i] with the index of the generator responsible for a cut,
4886  and place cuts flagged as global in the global cut pool for the model.
4887
4888  lastNumberCuts is the sum of cuts added in previous iterations; it's the
4889  offset to the proper starting position in whichGenerator.
4890*/
4891        int numberBefore =
4892          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
4893        int numberAfter =
4894          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
4895        // possibly extend whichGenerator
4896        resizeWhichGenerator(numberBefore, numberAfter);
4897        if (fullScan) {
4898          // counts
4899          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
4900        }
4901        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
4902       
4903        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
4904          whichGenerator_[numberBefore++] = i ;
4905          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
4906          if (thisCut->lb()>thisCut->ub())
4907            violated=-2; // sub-problem is infeasible
4908          if (thisCut->globallyValid()) {
4909            // add to global list
4910            OsiRowCut newCut(*thisCut);
4911            newCut.setGloballyValid(2);
4912            newCut.mutableRow().setTestForDuplicateIndex(false);
4913            globalCuts_.insert(newCut) ;
4914          }
4915        }
4916        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
4917          whichGenerator_[numberBefore++] = i ;
4918          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
4919          if (thisCut->globallyValid()) {
4920            // add to global list
4921            OsiColCut newCut(*thisCut);
4922            newCut.setGloballyValid(2);
4923            globalCuts_.insert(newCut) ;
4924          }
4925        }
4926      }
4927      // Add in any violated saved cuts
4928      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
4929        int numberOld = theseCuts.sizeRowCuts();
4930        int numberCuts = slackCuts.sizeRowCuts() ;
4931        int i;
4932        // possibly extend whichGenerator
4933        resizeWhichGenerator(numberOld, numberOld+numberCuts);
4934        for ( i = 0;i<numberCuts;i++) {
4935          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
4936          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
4937            if (messageHandler()->logLevel()>2)
4938              printf("Old cut added - violation %g\n",
4939                     thisCut->violated(cbcColSolution_)) ;
4940            whichGenerator_[numberOld++]=-1;
4941            theseCuts.insert(*thisCut) ;
4942          }
4943        }
4944      }
4945      delete [] eachCuts;
4946    }
4947    //if (!feasible)
4948    //break;
4949/*
4950  End of the loop to exercise each generator - try heuristics
4951  - unless at root node and first pass
4952*/
4953    if (numberNodes_||currentPassNumber_!=1) {
4954      double * newSolution = new double [numberColumns] ;
4955      double heuristicValue = getCutoff() ;
4956      int found = -1; // no solution found
4957      for (i = 0;i<numberHeuristics_;i++) {
4958        // see if heuristic will do anything
4959        double saveValue = heuristicValue ;
4960        int ifSol = 
4961          heuristic_[i]->solution(heuristicValue,
4962                                  newSolution,
4963                                  theseCuts) ;
4964        if (ifSol>0) {
4965          // better solution found
4966          found = i ;
4967          incrementUsed(newSolution);
4968        } else if (ifSol<0) {
4969          heuristicValue = saveValue ;
4970        }
4971      }
4972/*
4973  Did any of the heuristics turn up a new solution? Record it before we free
4974  the vector.
4975*/
4976      if (found >= 0) { 
4977        phase_=4;
4978        incrementUsed(newSolution);
4979        lastHeuristic_ = heuristic_[found];
4980        setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
4981        CbcTreeLocal * tree
4982          = dynamic_cast<CbcTreeLocal *> (tree_);
4983        if (tree)
4984          tree->passInSolution(bestSolution_,heuristicValue);
4985      }
4986      delete [] newSolution ;
4987    }
4988
4989#if 0
4990    // switch on to get all cuts printed
4991    theseCuts.printCuts() ;
4992#endif
4993    int numberColumnCuts = theseCuts.sizeColCuts() ;
4994    int numberRowCuts = theseCuts.sizeRowCuts() ;
4995    if (violated>=0)
4996      violated = numberRowCuts + numberColumnCuts ;
4997/*
4998  Apply column cuts (aka bound tightening). This may be partially redundant
4999  for column cuts returned by CglProbing, as generateCuts installs bounds
5000  from CglProbing when it determines it can fix a variable.
5001
5002  TODO: Looks like the use of violated has evolved. The value set above is
5003        completely ignored. All that's left is violated == -1 indicates some
5004        cut is violated, violated == -2 indicates infeasibility. Only
5005        infeasibility warrants exceptional action.
5006
5007  TODO: Strikes me that this code will fail to detect infeasibility, because
5008        the breaks escape the inner loops but the outer loop keeps going.
5009        Infeasibility in an early cut will be overwritten if a later cut is
5010        merely violated.
5011*/
5012    if (numberColumnCuts) {
5013
5014#ifdef CBC_DEBUG
5015      double * oldLower = new double [numberColumns] ;
5016      double * oldUpper = new double [numberColumns] ;
5017      memcpy(oldLower,cbcColLower_,numberColumns*sizeof(double)) ;
5018      memcpy(oldUpper,cbcColUpper_,numberColumns*sizeof(double)) ;
5019#endif
5020
5021      double integerTolerance = getDblParam(CbcIntegerTolerance) ;
5022      for (int i = 0;i<numberColumnCuts;i++) {
5023        const OsiColCut * thisCut = theseCuts.colCutPtr(i) ;
5024        const CoinPackedVector & lbs = thisCut->lbs() ;
5025        const CoinPackedVector & ubs = thisCut->ubs() ;
5026        int j ;
5027        int n ;
5028        const int * which ;
5029        const double * values ;
5030        n = lbs.getNumElements() ;
5031        which = lbs.getIndices() ;
5032        values = lbs.getElements() ;
5033        for (j = 0;j<n;j++) {
5034          int iColumn = which[j] ;
5035          double value = cbcColSolution_[iColumn] ;
5036#if CBC_DEBUG>1
5037          printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
5038                 cbcColSolution_[iColumn],oldUpper[iColumn],values[j]) ;
5039#endif
5040          solver_->setColLower(iColumn,values[j]) ;
5041          if (value<values[j]-integerTolerance)
5042            violated = -1 ;
5043          if (values[j]>cbcColUpper_[iColumn]+integerTolerance) {
5044            // infeasible
5045            violated = -2 ;
5046            break ;
5047          }
5048        }
5049        n = ubs.getNumElements() ;
5050        which = ubs.getIndices() ;
5051        values = ubs.getElements() ;
5052        for (j = 0;j<n;j++) {
5053          int iColumn = which[j] ;
5054          double value = cbcColSolution_[iColumn] ;
5055#if CBC_DEBUG>1
5056          printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
5057                 cbcColSolution_[iColumn],oldUpper[iColumn],values[j]) ;
5058#endif
5059          solver_->setColUpper(iColumn,values[j]) ;
5060          if (value>values[j]+integerTolerance)
5061            violated = -1 ;
5062          if (values[j]<cbcColLower_[iColumn]-integerTolerance) {
5063            // infeasible
5064            violated = -2 ;
5065            break ;
5066          }
5067        }
5068      }
5069#ifdef CBC_DEBUG
5070      delete [] oldLower ;
5071      delete [] oldUpper ;
5072#endif
5073    }
5074/*
5075  End installation of column cuts. The break here escapes the numberTries
5076  loop.
5077*/
5078    if (violated == -2||!feasible) {
5079      // infeasible
5080      feasible = false ;
5081      violated = -2;
5082      if (!numberNodes_) 
5083        messageHandler()->message(CBC_INFEAS,
5084                                  messages())
5085                                    << CoinMessageEol ;
5086      break ;
5087    }
5088/*
5089  Now apply the row (constraint) cuts. This is a bit more work because we need
5090  to obtain and augment the current basis.
5091
5092  TODO: Why do this work, if there are no row cuts? The current basis will do
5093        just fine.
5094*/
5095    int numberRowsNow = solver_->getNumRows() ;
5096    assert(numberRowsNow == numberRowsAtStart+lastNumberCuts) ;
5097    int numberToAdd = theseCuts.sizeRowCuts() ;
5098    numberNewCuts_ = lastNumberCuts+numberToAdd ;
5099/*
5100  Now actually add the row cuts and reoptimise.
5101
5102  Install the cuts in the solver using applyRowCuts and
5103  augment the basis with the corresponding slack. We also add each row cut to
5104  the set of row cuts (cuts.insert()) supplied as a parameter. The new basis
5105  must be set with setWarmStart().
5106
5107  TODO: Seems to me the original code could allocate addCuts with size 0, if
5108        numberRowCuts was 0 and numberColumnCuts was nonzero. That might
5109        explain the memory fault noted in the comment by AJK.  Unfortunately,
5110        just commenting out the delete[] results in massive memory leaks. Try
5111        a revision to separate the row cut case. Why do we need addCuts at
5112        all? A typing issue, apparently: OsiCut vs. OsiRowCut.
5113
5114  TODO: It looks to me as if numberToAdd and numberRowCuts are identical at
5115        this point. Confirm & get rid of one of them.
5116
5117  TODO: Any reason why the three loops can't be consolidated?
5118*/
5119    if (numberRowCuts > 0 || numberColumnCuts > 0)
5120    { if (numberToAdd > 0)
5121      { int i ;
5122        // Faster to add all at once
5123        const OsiRowCut ** addCuts = new const OsiRowCut * [numberToAdd] ;
5124        for (i = 0 ; i < numberToAdd ; i++)
5125        { addCuts[i] = &theseCuts.rowCut(i) ; }
5126        solver_->applyRowCuts(numberToAdd,addCuts) ;
5127        // AJK this caused a memory fault on Win32
5128        // may be okay now with ** form
5129        delete [] addCuts ;
5130        for (i = 0 ; i < numberToAdd ; i++)
5131        { cuts.insert(theseCuts.rowCut(i)) ; }
5132        CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5133        assert(basis != NULL); // make sure not volume
5134/* dylp bug
5135
5136  Consistent size used by OsiDylp as sanity check. Implicit resize seen
5137  as an error. Hence this call to resize is necessary.
5138*/
5139        basis->resize(numberRowsAtStart+numberNewCuts_,numberColumns) ;
5140        for (i = 0 ; i < numberToAdd ; i++) 
5141        { basis->setArtifStatus(numberRowsNow+i,
5142                                 CoinWarmStartBasis::basic) ; }
5143        if (solver_->setWarmStart(basis) == false)
5144        { throw CoinError("Fail setWarmStart() after cut installation.",
5145                          "solveWithCuts","CbcModel") ; }
5146        delete basis;
5147      }
5148      feasible = resolve(node ? node->nodeInfo() : NULL,2) ;
5149      if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
5150        numberTries=0; // exit
5151#     ifdef CBC_DEBUG
5152      printf("Obj value after cuts %g %d rows\n",solver_->getObjValue(),
5153              solver_->getNumRows()) ;
5154      if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
5155        assert(feasible) ;
5156#     endif
5157    }
5158/*
5159  No cuts. Cut short the cut generation (numberTries) loop.
5160*/
5161    else
5162    { numberTries = 0 ; }
5163/*
5164  If the problem is still feasible, first, call takeOffCuts() to remove cuts
5165  that are now slack. takeOffCuts() will call the solver to reoptimise if
5166  that's needed to restore a valid solution.
5167
5168  Next, see if we should quit due to diminishing returns:
5169    * we've tried three rounds of cut generation and we're getting
5170      insufficient improvement in the objective; or
5171    * we generated no cuts; or
5172    * the solver declared optimality with 0 iterations after we added the
5173      cuts generated in this round.
5174  If we decide to keep going, prep for the next iteration.
5175
5176  It just seems more safe to tell takeOffCuts() to call resolve(), even if
5177  we're not continuing cut generation. Otherwise code executed between here
5178  and final disposition of the node will need to be careful not to access the
5179  lp solution. It can happen that we lose feasibility in takeOffCuts ---
5180  numerical jitters when the cutoff bound is epsilon less than the current
5181  best, and we're evaluating an alternative optimum.
5182
5183  TODO: After successive rounds of code motion, there seems no need to
5184        distinguish between the two checks for aborting the cut generation
5185        loop. Confirm and clean up.
5186*/
5187    if (feasible)
5188    { int cutIterations = solver_->getIterationCount() ;
5189      if (numberOldActiveCuts_+numberNewCuts_) {
5190        OsiCuts * saveCuts = node ? NULL : &slackCuts;
5191        takeOffCuts(cuts,resolveAfterTakeOffCuts_,saveCuts) ;
5192        if (solver_->isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_)
5193          { feasible = false ;
5194#       ifdef CBC_DEBUG
5195          double z = solver_->getObjValue() ;
5196          double cut = getCutoff() ;
5197          printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n",
5198                 z-cut,z,cut) ;
5199#       endif
5200          }
5201      }
5202      if (feasible) {
5203        numberRowsAtStart = numberOldActiveCuts_+numberRowsAtContinuous_ ;
5204        lastNumberCuts = numberNewCuts_ ;
5205        if (direction*solver_->getObjValue() < lastObjective+minimumDrop &&
5206            currentPassNumber_ >= 3 && !keepGoing)
5207          { numberTries = 0 ; }
5208        if (numberRowCuts+numberColumnCuts == 0 || cutIterations == 0)
5209          { break ; }
5210        if (numberTries > 0) {
5211          reducedCostFix() ;
5212          lastObjective = direction*solver_->getObjValue() ;
5213        }
5214      }
5215    }
5216/*
5217  We've lost feasibility --- this node won't be referencing the cuts we've
5218  been collecting, so decrement the reference counts.
5219*/
5220    if (!feasible)
5221    { int i ;
5222      if (currentNumberCuts_) {
5223        lockThread();
5224        for (i = 0;i<currentNumberCuts_;i++) {
5225          // take off node
5226          if (addedCuts_[i]) {
5227            if (!addedCuts_[i]->decrement())
5228              delete addedCuts_[i] ;
5229            addedCuts_[i] = NULL ;
5230          }
5231        }
5232      }
5233      numberTries = 0 ;
5234    }
5235  } while (numberTries>0||keepGoing) ;
5236  {
5237    // switch on
5238    for (int i = 0;i<numberCutGenerators_;i++) 
5239      generator_[i]->setSwitchedOff(false);
5240  }
5241  //check feasibility.
5242  //If solution seems to be integer feasible calling setBestSolution
5243  //will eventually add extra global cuts which we need to install at
5244  //the nodes
5245
5246
5247  if(feasible&&solverCharacteristics_->solutionAddsCuts())  //check integer feasibility
5248  {
5249    bool integerFeasible = true;
5250    const double * save = testSolution_;
5251    testSolution_ = solver_->getColSolution();
5252    // point to useful information
5253    OsiBranchingInformation usefulInfo=usefulInformation();
5254    for (int i=0;i<numberObjects_ && integerFeasible;i++)
5255    {
5256      int preferredWay;
5257      double infeasibility = object_[i]->infeasibility(&usefulInfo,preferredWay);
5258      if(infeasibility)
5259        integerFeasible = false;
5260    }
5261    testSolution_ = save;
5262    if(integerFeasible)//update
5263    {
5264      double objValue = solver_->getObjValue();
5265      int numberGlobalBefore = globalCuts_.sizeRowCuts();
5266      // SOLUTION2 so won't up cutoff or print message
5267      setBestSolution(CBC_SOLUTION2, objValue, 
5268                      solver_->getColSolution(),0);
5269      int numberGlobalAfter = globalCuts_.sizeRowCuts();
5270      int numberToAdd = numberGlobalAfter - numberGlobalBefore;
5271      if(numberToAdd > 0)
5272        //We have added some cuts say they are tight at that node
5273        //Basis and lp should already have been updated
5274      {
5275        feasible = (solver_->isProvenOptimal() &&
5276                    !solver_->isDualObjectiveLimitReached()) ;
5277        if(feasible)
5278        {
5279          int numberCuts = numberNewCuts_ = cuts.sizeRowCuts();
5280      // possibly extend whichGenerator
5281          resizeWhichGenerator(numberCuts, numberToAdd+numberCuts);
5282         
5283          for (int i = numberGlobalBefore ; i < numberGlobalAfter ; i++)
5284            {       
5285              whichGenerator_[numberNewCuts_++]=-1;
5286#ifndef GLOBAL_CUTS_JUST_POINTERS
5287              cuts.insert(globalCuts_.rowCut(i)) ; 
5288#else
5289              OsiRowCut * rowCutPointer = globalCuts_.rowCutPtr(i);
5290              cuts.insert(rowCutPointer) ; 
5291#endif
5292            }
5293          numberNewCuts_ = lastNumberCuts+numberToAdd;
5294          //now take off the cuts which are not tight anymore
5295          takeOffCuts(cuts,resolveAfterTakeOffCuts_, NULL) ;
5296          if (solver_->isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_)
5297            { feasible = false ;}
5298        }
5299        if(!feasible) //node will be fathomed
5300          { 
5301          for (int i = 0;i<currentNumberCuts_;i++) 
5302            {
5303            // take off node
5304            if (addedCuts_[i])
5305              {
5306                if (!addedCuts_[i]->decrement())
5307                  delete addedCuts_[i] ;
5308                addedCuts_[i] = NULL ;
5309              }
5310            }
5311        }
5312      }
5313    }
5314  }
5315/*
5316  Reduced cost fix at end. Must also check feasible, in case we've popped out
5317  because a generator indicated we're infeasible.
5318*/
5319  if (feasible && solver_->isProvenOptimal())
5320    reducedCostFix();
5321  // If at root node do heuristics
5322  if (!numberNodes_) {
5323    // First see if any cuts are slack
5324    int numberRows = solver_->getNumRows();
5325    int numberAdded = numberRows-numberRowsAtContinuous_;
5326    if (numberAdded) {
5327      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5328      assert(basis != NULL);
5329      int * added = new int[numberAdded];
5330      int nDelete = 0;
5331      for (int j=numberRowsAtContinuous_;j<numberRows;j++) {
5332        if (basis->getArtifStatus(j)==CoinWarmStartBasis::basic) {
5333          //printf("%d slack!\n",j);
5334          added[nDelete++]=j;
5335        }
5336      }
5337      if (nDelete) {
5338        solver_->deleteRows(nDelete,added);
5339      }
5340      delete [] added;
5341      delete basis ;
5342    }
5343    // mark so heuristics can tell
5344    int savePass=currentPassNumber_;
5345    currentPassNumber_=999999;
5346    double * newSolution = new double [numberColumns] ;
5347    double heuristicValue = getCutoff() ;
5348    int found = -1; // no solution found
5349    for (int i = 0;i<numberHeuristics_;i++) {
5350      // see if heuristic will do anything
5351      double saveValue = heuristicValue ;
5352      int ifSol = heuristic_[i]->solution(heuristicValue,
5353                                          newSolution);
5354      if (ifSol>0) {
5355        // better solution found
5356        found = i ;
5357        incrementUsed(newSolution);
5358      } else {
5359        heuristicValue = saveValue ;
5360      }
5361    }
5362    currentPassNumber_=savePass;
5363    if (found >= 0) { 
5364      phase_=4;
5365      incrementUsed(newSolution);
5366      lastHeuristic_ = heuristic_[found];
5367      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
5368    }
5369    delete [] newSolution ;
5370  }
5371  // Up change due to cuts
5372  if (feasible)
5373    sumChangeObjective2_ += solver_->getObjValue()*solver_->getObjSense()
5374      - objectiveValue ;
5375  //if ((numberNodes_%100)==0)
5376  //printf("XXb sum obj changed by %g\n",sumChangeObjective2_);
5377/*
5378  End of cut generation loop.
5379
5380  Now, consider if we want to disable or adjust the frequency of use for any
5381  of the cut generators. If the client specified a positive number for
5382  howOften, it will never change. If the original value was negative, it'll
5383  be converted to 1000000+|howOften|, and this value will be adjusted each
5384  time fullScan is true. Actual cut generation is performed every
5385  howOften%1000000 nodes; the 1000000 offset is just a convenient way to
5386  specify that the frequency is adjustable.
5387
5388  During cut generation, we recorded the number of cuts produced by each
5389  generator for this node. For all cuts, whichGenerator records the generator
5390  that produced a cut.
5391
5392  TODO: All this should probably be hidden in a method of the CbcCutGenerator
5393  class.
5394
5395  TODO: Can the loop that scans over whichGenerator to accumulate per generator
5396        counts be replaced by values in countRowCuts and countColumnCuts?
5397*/
5398#ifdef NODE_LOG
5399  int fatherNum = (node == NULL) ? -1 : node->nodeNumber();
5400  double value =  (node == NULL) ? -1 : node->branchingObject()->value();
5401  string bigOne = (solver_->getIterationCount() > 30) ? "*******" : "";
5402  string way = (node == NULL) ? "" : (node->branchingObject()->way())==1 ? "Down" : "Up";
5403  std::cout<<"Node "<<numberNodes_<<", father "<<fatherNum<<", #iterations "<<solver_->getIterationCount()<<", sol value : "<<solver_->getObjValue()<<std::endl;
5404#endif
5405  if (fullScan&&numberCutGenerators_) {
5406    /* If cuts just at root node then it will probably be faster to
5407       update matrix and leave all in */
5408    int willBeCutsInTree=0;
5409    double thisObjective = solver_->getObjValue()*direction ;
5410    // get sizes
5411    int numberRowsAdded = solver_->getNumRows()-numberRowsAtStart;
5412    CoinBigIndex numberElementsAdded =  solver_->getNumElements()-numberElementsAtStart ;
5413    double densityOld = ((double) numberElementsAtStart)/((double) numberRowsAtStart);
5414    double densityNew = numberRowsAdded ? ((double) (numberElementsAdded))/((double) numberRowsAdded)
5415      : 0.0;
5416    if (!numberNodes_) {
5417      if (numberRowsAdded)
5418        handler_->message(CBC_CUTS_STATS,messages_)
5419          <<numberRowsAdded
5420          <<densityNew
5421          <<CoinMessageEol ;
5422      if (thisObjective-startObjective<1.0e-5&&numberElementsAdded>0.2*numberElementsAtStart)
5423        willBeCutsInTree=-1;
5424    }
5425    if ((numberRowsAdded>100+0.5*numberRowsAtStart
5426         ||numberElementsAdded>0.5*numberElementsAtStart)
5427        &&(densityNew>200.0&&numberRowsAdded>100&&densityNew>2.0*densityOld)) {
5428      // much bigger
5429      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
5430      //willBeCutsInTree=-1;
5431      //printf("Cuts will be taken off , %d rows added with density %g\n",
5432      //     numberRowsAdded,densityNew);
5433    }
5434    if (densityNew>100.0&&numberRowsAdded>2&&densityNew>2.0*densityOld) {
5435      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
5436      //willBeCutsInTree=-2;
5437      //printf("Density says no cuts ? , %d rows added with density %g\n",
5438      //     numberRowsAdded,densityNew);
5439    }
5440    // Root node or every so often - see what to turn off
5441    int i ;
5442    for (i = 0;i<numberCutGenerators_;i++) {
5443      int howOften = generator_[i]->howOften() ;
5444      if (howOften>-90) 
5445        willBeCutsInTree=0;
5446    }
5447    if (!numberNodes_)
5448      handler_->message(CBC_ROOT,messages_)
5449        <<numberNewCuts_
5450        <<startObjective<<thisObjective
5451        <<currentPassNumber_
5452        <<CoinMessageEol ;
5453    int * count = new int[numberCutGenerators_] ;
5454    memset(count,0,numberCutGenerators_*sizeof(int)) ;
5455    int numberActiveGenerators=0;
5456    for (i = 0;i<numberNewCuts_;i++) {
5457      int iGenerator = whichGenerator_[i];
5458      if (iGenerator>=0&&iGenerator<numberCutGenerators_)
5459        count[iGenerator]++ ;
5460    }
5461    double totalCuts = 0.0 ;
5462    //#define JUST_ACTIVE
5463    for (i = 0;i<numberCutGenerators_;i++) {
5464      if (countRowCuts[i]||countColumnCuts[i])
5465        numberActiveGenerators++;
5466#ifdef JUST_ACTIVE
5467      totalCuts += count[i] + 5.0*countColumnCuts[i] ;
5468#else
5469      totalCuts += countRowCuts[i] + 5.0*countColumnCuts[i] ;
5470#endif
5471    }
5472    int iProbing=-1;
5473    double small = (0.5* totalCuts) /
5474      ((double) numberActiveGenerators) ;
5475    for (i = 0;i<numberCutGenerators_;i++) {
5476      int howOften = generator_[i]->howOften() ;
5477      /*  Probing can be set to just do column cuts in treee.
5478          But if doing good then leave as on */
5479      bool probingWasOnBut=false;
5480      CglProbing * probing = dynamic_cast<CglProbing*>(generator_[i]->generator());
5481      if (probing&&!numberNodes_) {
5482        iProbing = i;
5483        if (probing->rowCuts()==-3) {
5484          probingWasOnBut=true;
5485          howOften = -98;
5486          probing->setRowCuts(3);
5487        }
5488      }
5489      if (willBeCutsInTree<0&&howOften==-98)
5490        howOften =-99;
5491      if (howOften==-98&&generator_[i]->switchOffIfLessThan()<0) {
5492        if (thisObjective-startObjective<0.005*fabs(startObjective)+1.0e-5)
5493          howOften=-99; // switch off
5494        if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5
5495            &&5*solver_->getNumRows()<solver_->getNumCols())
5496          howOften=-99; // switch off
5497      }
5498      if (howOften<-99)
5499        continue ;
5500      if (howOften<0||howOften >= 1000000) {
5501        if( !numberNodes_) {
5502          // If small number switch mostly off
5503#ifdef JUST_ACTIVE
5504          double thisCuts = count[i] + 5.0*countColumnCuts[i] ;
5505#else
5506          double thisCuts = countRowCuts[i] + 5.0*countColumnCuts[i] ;
5507#endif
5508          if (!thisCuts||howOften == -99) {
5509            if (howOften == -99||howOften == -98) 
5510              howOften = -100 ;
5511            else
5512              howOften = 1000000+SCANCUTS; // wait until next time
5513          } else if (thisCuts<small) {
5514            if (howOften!=1&&!probingWasOnBut) {
5515              if (generator_[i]->whatDepth()<0||howOften!=-1) {
5516                int k = (int) sqrt(small/thisCuts) ;
5517                if (howOften!=-98)
5518                  howOften = k+1000000 ;
5519                else
5520                  howOften=-100;
5521              } else {
5522                howOften=1;
5523              }
5524            } else {
5525              howOften=1;
5526              // allow cuts
5527              probingWasOnBut=false;
5528            }
5529          } else {
5530            if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5&&generator_[i]->whatDepth()<0)
5531              generator_[i]->setWhatDepth(5);
5532            howOften = 1+1000000 ;
5533          }
5534        }
5535        // If cuts useless switch off
5536        if (numberNodes_>=100000&&sumChangeObjective1_>2.0e2*(sumChangeObjective2_+1.0e-12)) {
5537          howOften = 1000000+SCANCUTS; // wait until next time
5538          //printf("switch off cut %d due to lack of use\n",i);
5539        }
5540      }
5541      if (!numberNodes_) {
5542        if (probingWasOnBut&&howOften==-100) {
5543          probing->setRowCuts(-3);
5544          howOften=1;
5545        }
5546        if (howOften==1)
5547          generator_[i]->setWhatDepth(1);
5548       
5549        if (howOften>=0&&generator_[i]->generator()->mayGenerateRowCutsInTree())
5550          willBeCutsInTree=1;
5551      }
5552      generator_[i]->setHowOften(howOften) ;
5553      if (howOften>=1000000&&howOften<2000000&&0) {
5554        // Go to depth
5555        int bias=1;
5556        if (howOften==1+1000000)
5557          generator_[i]->setWhatDepth(bias+1);
5558        else if (howOften<=10+1000000)
5559          generator_[i]->setWhatDepth(bias+2);
5560        else
5561          generator_[i]->setWhatDepth(bias+1000);
5562      }
5563      int newFrequency = generator_[i]->howOften()%1000000 ;
5564      // increment cut counts
5565      generator_[i]->incrementNumberCutsInTotal(countRowCuts[i]);
5566      generator_[i]->incrementNumberCutsActive(count[i]);
5567      if (handler_->logLevel()>1||!numberNodes_) {
5568        handler_->message(CBC_GENERATOR,messages_)
5569          <<i
5570          <<generator_[i]->cutGeneratorName()
5571          <<countRowCuts[i]<<count[i]
5572          <<countColumnCuts[i];
5573        handler_->printing(!numberNodes_&&generator_[i]->timing())
5574          <<generator_[i]->timeInCutGenerator();
5575        handler_->message()
5576          <<newFrequency
5577          <<CoinMessageEol ;
5578      }
5579    } 
5580    delete [] count ;
5581    if( !numberNodes_) {
5582      // save statistics
5583      for (i = 0;i<numberCutGenerators_;i++) {
5584        generator_[i]->setNumberCutsAtRoot(generator_[i]->numberCutsInTotal());
5585        generator_[i]->setNumberActiveCutsAtRoot(generator_[i]->numberCutsActive());
5586      }
5587      // decide on pseudo cost strategy
5588      int howOften = iProbing>=0 ? generator_[iProbing]->howOften() : 0;
5589      if ((howOften %1000000)!=1) 
5590        howOften = 0;
5591      //if (howOften) {
5592      //CglProbing * probing = dynamic_cast<CglProbing*>(generator_[iProbing]->generator());
5593      //}
5594      howOften = 0;
5595      if (howOften) {
5596        printf("** method 1\n");
5597        //CglProbing * probing = dynamic_cast<CglProbing*>(generator_[iProbing]->generator());
5598        generator_[iProbing]->setWhatDepth(1);
5599        // could set no row cuts
5600        //if (thisObjective-startObjective<0.001*fabs(startObjective)+1.0e-5)
5601        // probing->setRowCuts(0);
5602        for (int i=0;i<numberObjects_;i++) {
5603          CbcSimpleIntegerDynamicPseudoCost * obj =
5604            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
5605          if (obj) 
5606            obj->setMethod(1);
5607        }
5608      }
5609      if (willBeCutsInTree==-2)
5610        willBeCutsInTree=0;
5611      if( willBeCutsInTree<=0) {
5612        // Take off cuts
5613        cuts = OsiCuts();
5614        numberNewCuts_=0;
5615        if (!willBeCutsInTree) {
5616          // update size of problem
5617          numberRowsAtContinuous_ = solver_->getNumRows() ;
5618        } else {
5619          // take off cuts
5620          int numberRows = solver_->getNumRows();
5621          int numberAdded = numberRows-numberRowsAtContinuous_;
5622          if (numberAdded) {
5623            int * added = new int[numberAdded];
5624            for (int i=0;i<numberAdded;i++)
5625              added[i]=i+numberRowsAtContinuous_;
5626            solver_->deleteRows(numberAdded,added);
5627            delete [] added;
5628            // resolve so optimal
5629            resolve(solver_);
5630          }
5631        }
5632#ifdef COIN_HAS_CLP
5633        OsiClpSolverInterface * clpSolver
5634          = dynamic_cast<OsiClpSolverInterface *> (solver_);
5635        if (clpSolver) {
5636          // Maybe solver might like to know only column bounds will change
5637          //int options = clpSolver->specialOptions();
5638          //clpSolver->setSpecialOptions(options|128);
5639          clpSolver->synchronizeModel();
5640        }
5641#endif
5642      } else {
5643#ifdef COIN_HAS_CLP
5644        OsiClpSolverInterface * clpSolver
5645          = dynamic_cast<OsiClpSolverInterface *> (solver_);
5646        if (clpSolver) {
5647        // make sure factorization can't carry over
5648          int options = clpSolver->specialOptions();
5649          clpSolver->setSpecialOptions(options&(~8));
5650        }
5651#endif
5652      }
5653    }
5654  } else {
5655#ifdef COIN_HAS_CLP
5656    OsiClpSolverInterface * clpSolver
5657      = dynamic_cast<OsiClpSolverInterface *> (solver_);
5658    if (clpSolver) {
5659      // Maybe solver might like to know only column bounds will change
5660      //int options = clpSolver->specialOptions();
5661      //clpSolver->setSpecialOptions(options|128);
5662      clpSolver->synchronizeModel();
5663    }
5664#endif
5665    if (numberCutGenerators_) {
5666      int i;
5667      // add to counts anyway
5668      for (i = 0;i<numberCutGenerators_;i++) 
5669        generator_[i]->incrementNumberCutsInTotal(countRowCuts[i]);
5670      // What if not feasible as cuts may have helped
5671      if (feasible) {
5672        for (i = 0;i<numberNewCuts_;i++) {
5673          int iGenerator = whichGenerator_[i];
5674          if (iGenerator>=0)
5675            generator_[iGenerator]->incrementNumberCutsActive();
5676        }
5677      }
5678    }
5679  }
5680
5681  delete [] countRowCuts ;
5682  delete [] countColumnCuts ;
5683
5684
5685#ifdef CHECK_CUT_COUNTS
5686  if (feasible)
5687  { 
5688    CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5689    printf("solveWithCuts: Number of rows at end (only active cuts) %d\n",
5690           numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ;
5691    basis->print() ; delete basis;}
5692#endif
5693#ifdef CBC_DEBUG
5694  if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
5695    assert(feasible) ;
5696#endif
5697# ifdef COIN_HAS_CLP
5698  if (clpSolver) 
5699    clpSolver->setSpecialOptions(saveClpOptions);
5700# endif
5701#ifdef CBC_THREAD
5702  if (threadModel) {
5703    // stop threads
5704    int i;
5705    for (i=0;i<numberThreads_;i++) {
5706      pthread_cond_signal(threadInfo[i].condition2); // unlock
5707      while (threadInfo[i].returnCode==0) {
5708        struct timespec moreTime;
5709        moreTime.tv_nsec = 10000;
5710        moreTime.tv_sec = 0;
5711        nanosleep(&moreTime,NULL);
5712      }
5713      threadModel[i]->numberThreads_=0; // say exit
5714      threadInfo[i].returnCode=0;
5715      pthread_cond_signal(threadInfo[i].condition2); // unlock
5716      pthread_join(threadId[i],NULL);
5717      pthread_mutex_destroy (threadInfo[i].mutex2);
5718      pthread_cond_destroy (threadInfo[i].condition2);
5719      threadModel[i]->generator_[0]=NULL;
5720      delete [] threadModel[i]->generator_;
5721      threadModel[i]->generator_=NULL;
5722      threadModel[i]->object_=NULL;
5723    }
5724    pthread_cond_destroy (&condition_main);
5725    pthread_mutex_destroy (&condition_mutex);
5726    // delete models and solvers
5727    for (i=0;i<numberThreads_;i++) {
5728      delete threadModel[i];
5729    }
5730    delete [] mutex2;
5731    delete [] condition2;
5732    delete [] threadId;
5733    delete [] threadInfo;
5734    delete [] threadModel;
5735    mutex_ = saveMutex;
5736  }
5737#endif
5738
5739  return feasible ; }
5740
5741
5742/*
5743  Remove slack cuts. We obtain a basis and scan it. Cuts with basic slacks
5744  are purged. If any cuts are purged, resolve() is called to restore the
5745  solution held in the solver.  If resolve() pivots, there's the possibility
5746  that a slack may be pivoted in (trust me :-), so the process iterates.
5747  Setting allowResolve to false will suppress reoptimisation (but see note
5748  below).
5749
5750  At the level of the solver's constraint system, loose cuts are really
5751  deleted.  There's an implicit assumption that deleteRows will also update
5752  the active basis in the solver.
5753
5754  At the level of nodes and models, it's more complicated.
5755
5756  New cuts exist only in the collection of cuts passed as a parameter. They
5757  are deleted from the collection and that's the end of them.
5758
5759  Older cuts have made it into addedCuts_. Two separate actions are needed.
5760  The reference count for the CbcCountRowCut object is decremented. If this
5761  count falls to 0, the node which owns the cut is located, the reference to
5762  the cut is removed, and then the cut object is destroyed (courtesy of the
5763  CbcCountRowCut destructor). We also need to set the addedCuts_ entry to
5764  NULL. This is important so that when it comes time to generate basis edits
5765  we can tell this cut was dropped from the basis during processing of the
5766  node.
5767
5768  NOTE: In general, it's necessary to call resolve() after purging slack
5769        cuts.  Deleting constraints constitutes a change in the problem, and
5770        an OSI is not required to maintain a valid solution when the problem
5771        is changed. But ... it's really useful to maintain the active basis,
5772        and the OSI is supposed to do that. (Yes, it's splitting hairs.) In
5773        some places, it's possible to know that the solution will never be
5774        consulted after this call, only the basis.  (E.g., this routine is
5775        called as a last act before generating info to place the node in the
5776        live set.) For such use, set allowResolve to false.
5777 
5778  TODO: No real harm would be done if we just ignored the rare occasion when
5779        the call to resolve() pivoted a slack back into the basis. It's a
5780        minor inefficiency, at worst. But it does break assertions which
5781        check that there are no loose cuts in the basis. It might be better
5782        to remove the assertions.
5783*/
5784
5785void
5786CbcModel::takeOffCuts (OsiCuts &newCuts,
5787                       bool allowResolve, OsiCuts * saveCuts)
5788
5789{ // int resolveIterations = 0 ;
5790  int firstOldCut = numberRowsAtContinuous_ ;
5791  int totalNumberCuts = numberNewCuts_+numberOldActiveCuts_ ;
5792  int *solverCutIndices = new int[totalNumberCuts] ;
5793  int *newCutIndices = new int[numberNewCuts_] ;
5794  const CoinWarmStartBasis* ws ;
5795  CoinWarmStartBasis::Status status ;
5796  bool needPurge = true ;
5797/*
5798  The outer loop allows repetition of purge in the event that reoptimisation
5799  changes the basis. To start an iteration, clear the deletion counts and grab
5800  the current basis.
5801*/
5802  while (needPurge)
5803  { int numberNewToDelete = 0 ;
5804    int numberOldToDelete = 0 ;
5805    int i ;
5806    ws = dynamic_cast<const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5807/*
5808  Scan the basis entries of the old cuts generated prior to this round of cut
5809  generation.  Loose cuts are `removed' by decrementing their reference count
5810  and setting the addedCuts_ entry to NULL. (If the reference count falls to
5811  0, they're really deleted.  See CbcModel and CbcCountRowCut doc'n for
5812  principles of cut handling.)
5813*/
5814    int oldCutIndex = 0 ;
5815    if (numberOldActiveCuts_) {
5816      lockThread();
5817      for (i = 0 ; i < numberOldActiveCuts_ ; i++)
5818        { status = ws->getArtifStatus(i+firstOldCut) ;
5819        while (!addedCuts_[oldCutIndex]) oldCutIndex++ ;
5820        assert(oldCutIndex < currentNumberCuts_) ;
5821        // always leave if from nextRowCut_
5822        if (status == CoinWarmStartBasis::basic&&
5823            addedCuts_[oldCutIndex]->effectiveness()!=COIN_DBL_MAX)
5824          { solverCutIndices[numberOldToDelete++] = i+firstOldCut ;
5825          if (saveCuts) {
5826