source: branches/devel/Cbc/src/CbcNode.cpp @ 427

Last change on this file since 427 was 427, checked in by forrest, 13 years ago

fixes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 138.2 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
14#define CBC_WEAK_STRONG
15#include <cassert>
16#include <cfloat>
17#define CUTS
18#include "OsiSolverInterface.hpp"
19#include "OsiAuxInfo.hpp"
20#include "OsiSolverBranch.hpp"
21#include "CoinWarmStartBasis.hpp"
22#include "CoinTime.hpp"
23#include "CbcModel.hpp"
24#include "CbcNode.hpp"
25#include "CbcStatistics.hpp"
26#include "CbcStrategy.hpp"
27#include "CbcBranchActual.hpp"
28#include "CbcBranchDynamic.hpp"
29#include "OsiRowCut.hpp"
30#include "OsiRowCutDebugger.hpp"
31#include "OsiCuts.hpp"
32#include "CbcCountRowCut.hpp"
33#include "CbcFeasibilityBase.hpp"
34#include "CbcMessage.hpp"
35#ifdef COIN_HAS_CLP
36#include "OsiClpSolverInterface.hpp"
37#include "ClpSimplexOther.hpp"
38#endif
39using namespace std;
40#include "CglCutGenerator.hpp"
41// Default Constructor
42CbcNodeInfo::CbcNodeInfo ()
43  :
44  numberPointingToThis_(0),
45  parent_(NULL),
46  owner_(NULL),
47  numberCuts_(0),
48  nodeNumber_(0),
49  cuts_(NULL),
50  numberRows_(0),
51  numberBranchesLeft_(0)
52{
53#ifdef CHECK_NODE
54  printf("CbcNodeInfo %x Constructor\n",this);
55#endif
56}
57// Constructor given parent
58CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent)
59  :
60  numberPointingToThis_(2),
61  parent_(parent),
62  owner_(NULL),
63  numberCuts_(0),
64  nodeNumber_(0),
65  cuts_(NULL),
66  numberRows_(0),
67  numberBranchesLeft_(2)
68{
69#ifdef CHECK_NODE
70  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
71#endif
72  if (parent_) {
73    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
74  }
75}
76// Copy Constructor
77CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs)
78  :
79  numberPointingToThis_(rhs.numberPointingToThis_),
80  parent_(rhs.parent_),
81  owner_(rhs.owner_),
82  numberCuts_(rhs.numberCuts_),
83  nodeNumber_(rhs.nodeNumber_),
84  cuts_(NULL),
85  numberRows_(rhs.numberRows_),
86  numberBranchesLeft_(rhs.numberBranchesLeft_)
87{
88#ifdef CHECK_NODE
89  printf("CbcNodeInfo %x Copy constructor\n",this);
90#endif
91  if (numberCuts_) {
92    cuts_ = new CbcCountRowCut * [numberCuts_];
93    int n=0;
94    for (int i=0;i<numberCuts_;i++) {
95      CbcCountRowCut * thisCut = rhs.cuts_[i];
96      if (thisCut) {
97        // I think this is correct - new one should take priority
98        thisCut->setInfo(this,n);
99        thisCut->increment(numberBranchesLeft_); 
100        cuts_[n++] = thisCut;
101      }
102    }
103    numberCuts_=n;
104  }
105}
106// Constructor given parent and owner
107CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent, CbcNode * owner)
108  :
109  numberPointingToThis_(2),
110  parent_(parent),
111  owner_(owner),
112  numberCuts_(0),
113  nodeNumber_(0),
114  cuts_(NULL),
115  numberRows_(0),
116  numberBranchesLeft_(2)
117{
118#ifdef CHECK_NODE
119  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
120#endif
121  if (parent_) {
122    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
123  }
124}
125
126/**
127  Take care to detach from the owning CbcNode and decrement the reference
128  count in the parent.  If this is the last nodeInfo object pointing to the
129  parent, make a recursive call to delete the parent.
130*/
131CbcNodeInfo::~CbcNodeInfo()
132{
133#ifdef CHECK_NODE
134  printf("CbcNodeInfo %x Destructor parent %x\n",this,parent_);
135#endif
136
137  assert(!numberPointingToThis_);
138  // But they may be some left (max nodes?)
139  for (int i=0;i<numberCuts_;i++) 
140    delete cuts_[i];
141
142  delete [] cuts_;
143  if (owner_) 
144    owner_->nullNodeInfo();
145  if (parent_) {
146    int numberLinks = parent_->decrement();
147    if (!numberLinks) delete parent_;
148  }
149}
150
151
152//#define ALLCUTS
153void
154CbcNodeInfo::decrementCuts(int change)
155{
156  int i;
157  // get rid of all remaining if negative
158  int changeThis;
159  if (change<0)
160    changeThis = numberBranchesLeft_;
161  else
162    changeThis = change;
163 // decrement cut counts
164  for (i=0;i<numberCuts_;i++) {
165    if (cuts_[i]) {
166      int number = cuts_[i]->decrement(changeThis);
167      if (!number) {
168        //printf("info %x del cut %d %x\n",this,i,cuts_[i]);
169        delete cuts_[i];
170        cuts_[i]=NULL;
171      }
172    }
173  }
174}
175void
176CbcNodeInfo::decrementParentCuts(int change)
177{
178  if (parent_) {
179    // get rid of all remaining if negative
180    int changeThis;
181    if (change<0)
182      changeThis = numberBranchesLeft_;
183    else
184      changeThis = change;
185    int i;
186    // Get over-estimate of space needed for basis
187    CoinWarmStartBasis dummy;
188    dummy.setSize(0,numberRows_+numberCuts_);
189    buildRowBasis(dummy);
190    /* everything is zero (i.e. free) so we can use to see
191       if latest basis */
192    CbcNodeInfo * thisInfo = parent_;
193    while (thisInfo) 
194      thisInfo = thisInfo->buildRowBasis(dummy);
195    // decrement cut counts
196    thisInfo = parent_;
197    int numberRows=numberRows_;
198    while (thisInfo) {
199      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
200        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
201#ifdef ALLCUTS
202        status = CoinWarmStartBasis::isFree;
203#endif
204        if (thisInfo->cuts_[i]) {
205          int number=1;
206          if (status!=CoinWarmStartBasis::basic) {
207            // tight - drop 1 or 2
208            if (change<0)
209              number = thisInfo->cuts_[i]->decrement(changeThis);
210            else
211              number = thisInfo->cuts_[i]->decrement(change);
212          }
213          if (!number) {
214            delete thisInfo->cuts_[i];
215            thisInfo->cuts_[i]=NULL;
216          }
217        }
218      }
219      thisInfo = thisInfo->parent_;
220    }
221  }
222}
223
224void
225CbcNodeInfo::incrementParentCuts(int change)
226{
227  if (parent_) {
228    int i;
229    // Get over-estimate of space needed for basis
230    CoinWarmStartBasis dummy;
231    dummy.setSize(0,numberRows_+numberCuts_);
232    /* everything is zero (i.e. free) so we can use to see
233       if latest basis */
234    buildRowBasis(dummy);
235    CbcNodeInfo * thisInfo = parent_;
236    while (thisInfo) 
237      thisInfo = thisInfo->buildRowBasis(dummy);
238    // increment cut counts
239    thisInfo = parent_;
240    int numberRows=numberRows_;
241    while (thisInfo) {
242      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
243        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
244#ifdef ALLCUTS
245        status = CoinWarmStartBasis::isFree;
246#endif
247        if (thisInfo->cuts_[i]&&status!=CoinWarmStartBasis::basic) {
248          thisInfo->cuts_[i]->increment(change);
249        }
250      }
251      thisInfo = thisInfo->parent_;
252    }
253  }
254}
255
256/*
257  Append cuts to the cuts_ array in a nodeInfo. The initial reference count
258  is set to numberToBranchOn, which will normally be the number of arms
259  defined for the CbcBranchingObject attached to the CbcNode that owns this
260  CbcNodeInfo.
261*/
262void
263CbcNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn,
264                      int * whichGenerator)
265{
266  int numberCuts = cuts.sizeRowCuts();
267  if (numberCuts) {
268    int i;
269    if (!numberCuts_) {
270      cuts_ = new CbcCountRowCut * [numberCuts];
271    } else {
272      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
273      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
274      delete [] cuts_;
275      cuts_ = temp;
276    }
277    for (i=0;i<numberCuts;i++) {
278      CbcCountRowCut * thisCut = new CbcCountRowCut(*cuts.rowCutPtr(i),
279                                                    this,numberCuts_);
280      thisCut->increment(numberToBranchOn); 
281      cuts_[numberCuts_++] = thisCut;
282#ifdef CBC_DEBUG
283#if CBC_DEBUG>1
284      int n=thisCut->row().getNumElements();
285      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
286             thisCut->ub());
287      int j;
288      const int * index = thisCut->row().getIndices();
289      const double * element = thisCut->row().getElements();
290      for (j=0;j<n;j++) {
291        printf(" (%d,%g)",index[j],element[j]);
292        assert(fabs(element[j])>1.00e-12);
293      }
294      printf("\n");
295#else
296      int n=thisCut->row().getNumElements();
297      int j;
298      const double * element = thisCut->row().getElements();
299      for (j=0;j<n;j++) {
300        assert(fabs(element[j])>1.00e-12);
301      }
302#endif
303#endif
304    }
305  }
306}
307
308void
309CbcNodeInfo::addCuts(int numberCuts, CbcCountRowCut ** cut, 
310                     int numberToBranchOn)
311{
312  if (numberCuts) {
313    int i;
314    if (!numberCuts_) {
315      cuts_ = new CbcCountRowCut * [numberCuts];
316    } else {
317      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
318      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
319      delete [] cuts_;
320      cuts_ = temp;
321    }
322    for (i=0;i<numberCuts;i++) {
323      CbcCountRowCut * thisCut = cut[i];
324      thisCut->setInfo(this,numberCuts_);
325      //printf("info %x cut %d %x\n",this,i,thisCut);
326      thisCut->increment(numberToBranchOn); 
327      cuts_[numberCuts_++] = thisCut;
328#ifdef CBC_DEBUG
329      int n=thisCut->row().getNumElements();
330#if CBC_DEBUG>1
331      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
332             thisCut->ub());
333#endif
334      int j;
335#if CBC_DEBUG>1
336      const int * index = thisCut->row().getIndices();
337#endif
338      const double * element = thisCut->row().getElements();
339      for (j=0;j<n;j++) {
340#if CBC_DEBUG>1
341        printf(" (%d,%g)",index[j],element[j]);
342#endif
343        assert(fabs(element[j])>1.00e-12);
344      }
345      printf("\n");
346#endif
347    }
348  }
349}
350
351// delete cuts
352void
353CbcNodeInfo::deleteCuts(int numberToDelete, CbcCountRowCut ** cuts)
354{
355  int i;
356  int j;
357  int last=-1;
358  for (i=0;i<numberToDelete;i++) {
359    CbcCountRowCut * next = cuts[i];
360    for (j=last+1;j<numberCuts_;j++) {
361      if (next==cuts_[j])
362        break;
363    }
364    if (j==numberCuts_) {
365      // start from beginning
366      for (j=0;j<last;j++) {
367        if (next==cuts_[j])
368          break;
369      }
370      assert(j<last);
371    }
372    last=j;
373    int number = cuts_[j]->decrement();
374    if (!number)
375      delete cuts_[j];
376    cuts_[j]=NULL;
377  }
378  j=0;
379  for (i=0;i<numberCuts_;i++) {
380    if (cuts_[i])
381      cuts_[j++]=cuts_[i];
382  }
383  numberCuts_ = j;
384}
385
386// delete cuts
387void
388CbcNodeInfo::deleteCuts(int numberToDelete, int * which)
389{
390  int i;
391  for (i=0;i<numberToDelete;i++) {
392    int iCut=which[i];
393    int number = cuts_[iCut]->decrement();
394    if (!number)
395      delete cuts_[iCut];
396    cuts_[iCut]=NULL;
397  }
398  int j=0;
399  for (i=0;i<numberCuts_;i++) {
400    if (cuts_[i])
401      cuts_[j++]=cuts_[i];
402  }
403  numberCuts_ = j;
404}
405
406// Really delete a cut
407void 
408CbcNodeInfo::deleteCut(int whichOne)
409{
410  assert(whichOne<numberCuts_);
411  cuts_[whichOne]=NULL;
412}
413
414CbcFullNodeInfo::CbcFullNodeInfo() :
415  CbcNodeInfo(),
416  basis_(),
417  numberIntegers_(0),
418  lower_(NULL),
419  upper_(NULL)
420{
421}
422CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model,
423                                 int numberRowsAtContinuous) :
424  CbcNodeInfo()
425{
426  OsiSolverInterface * solver = model->solver();
427  numberRows_ = numberRowsAtContinuous;
428  numberIntegers_ = model->numberIntegers();
429  int numberColumns = model->getNumCols();
430  lower_ = new double [numberColumns];
431  upper_ = new double [numberColumns];
432  const double * lower = solver->getColLower();
433  const double * upper = solver->getColUpper();
434  int i;
435
436  for (i=0;i<numberColumns;i++) {
437    lower_[i]=lower[i];
438    upper_[i]=upper[i];
439  }
440
441  basis_ =  dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
442}
443
444CbcFullNodeInfo::CbcFullNodeInfo(const CbcFullNodeInfo & rhs) :
445  CbcNodeInfo(rhs)
446{ 
447  basis_= dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
448  numberIntegers_=rhs.numberIntegers_;
449  lower_=NULL;
450  upper_=NULL;
451  if (rhs.lower_!=NULL) {
452    int numberColumns = basis_->getNumStructural();
453    lower_ = new double [numberColumns];
454    upper_ = new double [numberColumns];
455    assert (upper_!=NULL);
456    memcpy(lower_,rhs.lower_,numberColumns*sizeof(double));
457    memcpy(upper_,rhs.upper_,numberColumns*sizeof(double));
458  }
459}
460
461CbcNodeInfo * 
462CbcFullNodeInfo::clone() const
463{ 
464  return (new CbcFullNodeInfo(*this));
465}
466
467CbcFullNodeInfo::~CbcFullNodeInfo ()
468{
469  delete basis_ ;
470  delete [] lower_;
471  delete [] upper_;
472}
473
474/*
475   The basis supplied as a parameter is deleted and replaced with a new basis
476   appropriate for the node, and lower and upper bounds on variables are
477   reset according to the stored bounds arrays. Any cuts associated with this
478   node are added to the list in addCuts, but not actually added to the
479   constraint system in the model.
480
481   Why pass in a basis at all? The short answer is ``We need the parameter to
482   pass out a basis, so might as well use it to pass in the size.''
483   
484   A longer answer is that in practice we take a memory allocation hit up in
485   addCuts1 (the only place applyToModel is called) when we setSize() the
486   basis that's passed in. It's immediately tossed here in favour of a clone
487   of the basis attached to this nodeInfo. This can probably be fixed, given
488   a bit of thought.
489*/
490
491void CbcFullNodeInfo::applyToModel (CbcModel *model,
492                                    CoinWarmStartBasis *&basis,
493                                    CbcCountRowCut **addCuts,
494                                    int &currentNumberCuts) const 
495
496{ OsiSolverInterface *solver = model->solver() ;
497
498  // branch - do bounds
499  int i;
500  solver->setColLower(lower_);
501  solver->setColUpper(upper_);
502  int numberColumns = model->getNumCols();
503  // move basis - but make sure size stays
504  // for bon-min - should not be needed int numberRows = model->getNumRows();
505  int numberRows=basis->getNumArtificial();
506  delete basis ;
507  if (basis_) {
508    basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
509    basis->resize(numberRows,numberColumns);
510  } else {
511    // We have a solver without a basis
512    basis=NULL;
513  }
514  for (i=0;i<numberCuts_;i++) 
515    addCuts[currentNumberCuts+i]= cuts_[i];
516  currentNumberCuts += numberCuts_;
517  assert(!parent_);
518  return ;
519}
520
521/* Builds up row basis backwards (until original model).
522   Returns NULL or previous one to apply .
523   Depends on Free being 0 and impossible for cuts
524*/
525CbcNodeInfo * 
526CbcFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
527{
528  const unsigned int * saved = 
529    (const unsigned int *) basis_->getArtificialStatus();
530  unsigned int * now = 
531    (unsigned int *) basis.getArtificialStatus();
532  int number=basis_->getNumArtificial()>>4;;
533  int i;
534  for (i=0;i<number;i++) { 
535    if (!now[i])
536      now[i] = saved[i];
537  }
538  return NULL;
539}
540
541// Default constructor
542CbcPartialNodeInfo::CbcPartialNodeInfo()
543
544  : CbcNodeInfo(),
545    basisDiff_(NULL),
546    variables_(NULL),
547    newBounds_(NULL),
548    numberChangedBounds_(0)
549
550{ /* this space intentionally left blank */ }
551
552// Constructor from current state
553CbcPartialNodeInfo::CbcPartialNodeInfo (CbcNodeInfo *parent, CbcNode *owner,
554                                        int numberChangedBounds,
555                                        const int *variables,
556                                        const double *boundChanges,
557                                        const CoinWarmStartDiff *basisDiff)
558 : CbcNodeInfo(parent,owner)
559{
560  basisDiff_ = basisDiff->clone() ;
561
562  numberChangedBounds_ = numberChangedBounds;
563  variables_ = new int [numberChangedBounds_];
564  newBounds_ = new double [numberChangedBounds_];
565
566  int i ;
567  for (i=0;i<numberChangedBounds_;i++) {
568    variables_[i]=variables[i];
569    newBounds_[i]=boundChanges[i];
570  }
571}
572
573CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs)
574
575  : CbcNodeInfo(rhs.parent_)
576
577{ basisDiff_ = rhs.basisDiff_->clone() ;
578
579  numberChangedBounds_ = rhs.numberChangedBounds_;
580  variables_ = new int [numberChangedBounds_];
581  newBounds_ = new double [numberChangedBounds_];
582
583  int i ;
584  for (i=0;i<numberChangedBounds_;i++) {
585    variables_[i]=rhs.variables_[i];
586    newBounds_[i]=rhs.newBounds_[i];
587  }
588}
589
590CbcNodeInfo * 
591CbcPartialNodeInfo::clone() const
592{ 
593  return (new CbcPartialNodeInfo(*this));
594}
595
596
597CbcPartialNodeInfo::~CbcPartialNodeInfo ()
598{
599  delete basisDiff_ ;
600  delete [] variables_;
601  delete [] newBounds_;
602}
603
604
605/**
606   The basis supplied as a parameter is incrementally modified, and lower and
607   upper bounds on variables in the model are incrementally modified. Any
608   cuts associated with this node are added to the list in addCuts.
609*/
610
611void CbcPartialNodeInfo::applyToModel (CbcModel *model,
612                                       CoinWarmStartBasis *&basis,
613                                       CbcCountRowCut **addCuts,
614                                       int &currentNumberCuts) const 
615
616{ OsiSolverInterface *solver = model->solver();
617
618  basis->applyDiff(basisDiff_) ;
619
620  // branch - do bounds
621  int i;
622  for (i=0;i<numberChangedBounds_;i++) {
623    int variable = variables_[i];
624    if ((variable&0x80000000)==0) {
625      // lower bound changing
626      solver->setColLower(variable,newBounds_[i]);
627    } else {
628      // upper bound changing
629      solver->setColUpper(variable&0x7fffffff,newBounds_[i]);
630    }
631  }
632  for (i=0;i<numberCuts_;i++) {
633    addCuts[currentNumberCuts+i]= cuts_[i];
634    if (cuts_[i]&&model->messageHandler()->logLevel()>4) {
635      cuts_[i]->print();
636    }
637  }
638   
639  currentNumberCuts += numberCuts_;
640  return ;
641}
642
643/* Builds up row basis backwards (until original model).
644   Returns NULL or previous one to apply .
645   Depends on Free being 0 and impossible for cuts
646*/
647
648CbcNodeInfo * 
649CbcPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
650
651{ basis.applyDiff(basisDiff_) ;
652
653  return parent_ ; }
654
655
656CbcNode::CbcNode() :
657  nodeInfo_(NULL),
658  objectiveValue_(1.0e100),
659  guessedObjectiveValue_(1.0e100),
660  branch_(NULL),
661  depth_(-1),
662  numberUnsatisfied_(0)
663{
664#ifdef CHECK_NODE
665  printf("CbcNode %x Constructor\n",this);
666#endif
667}
668
669CbcNode::CbcNode(CbcModel * model,
670                 CbcNode * lastNode) :
671  nodeInfo_(NULL),
672  objectiveValue_(1.0e100),
673  guessedObjectiveValue_(1.0e100),
674  branch_(NULL),
675  depth_(-1),
676  numberUnsatisfied_(0)
677{
678#ifdef CHECK_NODE
679  printf("CbcNode %x Constructor from model\n",this);
680#endif
681  model->setObjectiveValue(this,lastNode);
682
683  if (lastNode)
684    if (lastNode->nodeInfo_)
685    lastNode->nodeInfo_->increment();
686}
687
688
689void
690CbcNode::createInfo (CbcModel *model,
691                     CbcNode *lastNode,
692                     const CoinWarmStartBasis *lastws,
693                     const double *lastLower, const double *lastUpper,
694                     int numberOldActiveCuts,int numberNewCuts)
695{ OsiSolverInterface * solver = model->solver();
696 CbcStrategy * strategy = model->strategy();
697/*
698  The root --- no parent. Create full basis and bounds information.
699*/
700  if (!lastNode)
701  { 
702    if (!strategy)
703      nodeInfo_=new CbcFullNodeInfo(model,solver->getNumRows());
704    else
705      nodeInfo_ = strategy->fullNodeInfo(model,solver->getNumRows());
706  }
707/*
708  Not the root. Create an edit from the parent's basis & bound information.
709  This is not quite as straightforward as it seems. We need to reintroduce
710  cuts we may have dropped out of the basis, in the correct position, because
711  this whole process is strictly positional. Start by grabbing the current
712  basis.
713*/
714  else
715  { const CoinWarmStartBasis* ws =
716      dynamic_cast<const CoinWarmStartBasis*>(solver->getWarmStart());
717    assert(ws!=NULL); // make sure not volume
718    //int numberArtificials = lastws->getNumArtificial();
719    int numberColumns = solver->getNumCols();
720   
721    const double * lower = solver->getColLower();
722    const double * upper = solver->getColUpper();
723
724    int i;
725/*
726  Create a clone and resize it to hold all the structural constraints, plus
727  all the cuts: old cuts, both active and inactive (currentNumberCuts), and
728  new cuts (numberNewCuts).
729
730  TODO: You'd think that the set of constraints (logicals) in the expanded
731        basis should match the set represented in lastws. At least, that's
732        what I thought. But at the point I first looked hard at this bit of
733        code, it turned out that lastws was the stripped basis produced at
734        the end of addCuts(), rather than the raw basis handed back by
735        addCuts1(). The expanded basis here is equivalent to the raw basis of
736        addCuts1(). I said ``whoa, that's not good, I must have introduced a
737        bug'' and went back to John's code to see where I'd gone wrong.
738        And discovered the same `error' in his code.
739
740        After a bit of thought, my conclusion is that correctness is not
741        affected by whether lastws is the stripped or raw basis. The diffs
742        have no semantics --- just a set of changes that need to be made
743        to convert lastws into expanded. I think the only effect is that we
744        store a lot more diffs (everything in expanded that's not covered by
745        the stripped basis). But I need to give this more thought. There
746        may well be some subtle error cases.
747
748        In the mean time, I've twiddled addCuts() to set lastws to the raw
749        basis. Makes me (Lou) less nervous to compare apples to apples.
750*/
751    CoinWarmStartBasis *expanded = 
752      dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
753    int numberRowsAtContinuous = model->numberRowsAtContinuous();
754    int iFull = numberRowsAtContinuous+model->currentNumberCuts()+
755      numberNewCuts;
756    //int numberArtificialsNow = iFull;
757    //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
758    //printf("l %d full %d\n",maxBasisLength,iFull);
759    if (expanded) 
760      expanded->resize(iFull,numberColumns);
761#ifdef CBC_CHECK_BASIS
762    printf("Before expansion: orig %d, old %d, new %d, current %d\n",
763           numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts,
764           model->currentNumberCuts()) ;
765    ws->print();
766#endif
767/*
768  Now fill in the expanded basis. Any indices beyond nPartial must
769  be cuts created while processing this node --- they can be copied directly
770  into the expanded basis. From nPartial down, pull the status of active cuts
771  from ws, interleaving with a B entry for the deactivated (loose) cuts.
772*/
773    int numberDropped = model->currentNumberCuts()-numberOldActiveCuts;
774    int iCompact=iFull-numberDropped;
775    CbcCountRowCut ** cut = model->addedCuts();
776    int nPartial = model->currentNumberCuts()+numberRowsAtContinuous;
777    iFull--;
778    for (;iFull>=nPartial;iFull--) {
779      CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
780      //assert (status != CoinWarmStartBasis::basic); // may be permanent cut
781      expanded->setArtifStatus(iFull,status);
782    }
783    for (;iFull>=numberRowsAtContinuous;iFull--) {
784      if (cut[iFull-numberRowsAtContinuous]) {
785        CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
786        // If no cut generator being used then we may have basic variables
787        //if (model->getMaximumCutPasses()&&
788        //  status == CoinWarmStartBasis::basic)
789        //printf("cut basic\n");
790        expanded->setArtifStatus(iFull,status);
791      } else {
792        expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic);
793      }
794    }
795#ifdef CBC_CHECK_BASIS
796    printf("Expanded basis\n");
797    expanded->print() ;
798    printf("Diffing against\n") ;
799    lastws->print() ;
800#endif   
801/*
802  Now that we have two bases in proper positional correspondence, creating
803  the actual diff is dead easy.
804*/
805
806    CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
807/*
808  Diff the bound vectors. It's assumed the number of structural variables is
809  not changing. Assuming that branching objects all involve integer variables,
810  we should see at least one bound change as a consequence of processing this
811  subproblem. Different types of branching objects could break this assertion.
812  Not true at all - we have not applied current branch - JJF.
813*/
814    double *boundChanges = new double [2*numberColumns] ;
815    int *variables = new int [2*numberColumns] ;
816    int numberChangedBounds=0;
817    for (i=0;i<numberColumns;i++) {
818      if (lower[i]!=lastLower[i]) {
819        variables[numberChangedBounds]=i;
820        boundChanges[numberChangedBounds++]=lower[i];
821      }
822      if (upper[i]!=lastUpper[i]) {
823        variables[numberChangedBounds]=i|0x80000000;
824        boundChanges[numberChangedBounds++]=upper[i];
825      }
826#ifdef CBC_DEBUG
827      if (lower[i]!=lastLower[i])
828        printf("lower on %d changed from %g to %g\n",
829               i,lastLower[i],lower[i]);
830      if (upper[i]!=lastUpper[i])
831        printf("upper on %d changed from %g to %g\n",
832               i,lastUpper[i],upper[i]);
833#endif
834    }
835#ifdef CBC_DEBUG
836    printf("%d changed bounds\n",numberChangedBounds) ;
837#endif
838    //if (lastNode->branchingObject()->boundBranch())
839    //assert (numberChangedBounds);
840/*
841  Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
842  return.
843*/
844    if (!strategy)
845      nodeInfo_ =
846        new CbcPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds,
847                               variables,boundChanges,basisDiff) ;
848    else
849      nodeInfo_ = strategy->partialNodeInfo(model, lastNode->nodeInfo_,this,numberChangedBounds,
850                               variables,boundChanges,basisDiff) ;
851    delete basisDiff ;
852    delete [] boundChanges;
853    delete [] variables;
854    delete expanded ;
855    delete ws;
856  }
857  // Set node number
858  nodeInfo_->setNodeNumber(model->getNodeCount2());
859}
860
861/*
862  The routine scans through the object list of the model looking for objects
863  that indicate infeasibility. It tests each object using strong branching
864  and selects the one with the least objective degradation.  A corresponding
865  branching object is left attached to lastNode.
866
867  If strong branching is disabled, a candidate object is chosen essentially
868  at random (whatever object ends up in pos'n 0 of the candidate array).
869
870  If a branching candidate is found to be monotone, bounds are set to fix the
871  variable and the routine immediately returns (the caller is expected to
872  reoptimize).
873
874  If a branching candidate is found to result in infeasibility in both
875  directions, the routine immediately returns an indication of infeasibility.
876
877  Returns:  0   both branch directions are feasible
878           -1   branching variable is monotone
879           -2   infeasible
880
881  Original comments:
882    Here could go cuts etc etc
883    For now just fix on objective from strong branching.
884*/
885
886int CbcNode::chooseBranch (CbcModel *model, CbcNode *lastNode,int numberPassesLeft)
887
888{ if (lastNode)
889    depth_ = lastNode->depth_+1;
890  else
891    depth_ = 0;
892  delete branch_;
893  branch_=NULL;
894  OsiSolverInterface * solver = model->solver();
895  double saveObjectiveValue = solver->getObjValue();
896  double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
897  const double * lower = solver->getColLower();
898  const double * upper = solver->getColUpper();
899  // See what user thinks
900  int anyAction=model->problemFeasibility()->feasible(model,0);
901  if (anyAction) {
902    // will return -2 if infeasible , 0 if treat as integer
903    return anyAction-1;
904  }
905  double integerTolerance = 
906    model->getDblParam(CbcModel::CbcIntegerTolerance);
907  int i;
908  bool beforeSolution = model->getSolutionCount()==0;
909  int numberStrong=model->numberStrong();
910  // switch off strong if hotstart
911  if (model->hotstartSolution())
912    numberStrong=0;
913  int numberStrongDone=0;
914  int numberUnfinished=0;
915  int numberStrongInfeasible=0;
916  int numberStrongIterations=0;
917  int saveNumberStrong=numberStrong;
918  int numberObjects = model->numberObjects();
919  bool checkFeasibility = numberObjects>model->numberIntegers();
920  int maximumStrong = CoinMax(CoinMin(model->numberStrong(),numberObjects),1);
921  int numberColumns = model->getNumCols();
922  double * saveUpper = new double[numberColumns];
923  double * saveLower = new double[numberColumns];
924
925  // Save solution in case heuristics need good solution later
926 
927  double * saveSolution = new double[numberColumns];
928  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
929  model->reserveCurrentSolution(saveSolution);
930  /*
931    Get a branching decision object. Use the default decision criteria unless
932    the user has loaded a decision method into the model.
933  */
934  CbcBranchDecision *decision = model->branchingMethod();
935  CbcDynamicPseudoCostBranchingObject * dynamicBranchingObject =
936    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(decision);
937  if (!decision||dynamicBranchingObject)
938    decision = new CbcBranchDefaultDecision();
939
940  CbcStrongInfo * choice = new CbcStrongInfo[maximumStrong];
941  for (i=0;i<numberColumns;i++) {
942    saveLower[i] = lower[i];
943    saveUpper[i] = upper[i];
944  }
945  // May go round twice if strong branching fixes all local candidates
946  bool finished=false;
947  double estimatedDegradation=0.0; 
948  while(!finished) {
949    finished=true;
950    // Some objects may compute an estimate of best solution from here
951    estimatedDegradation=0.0; 
952    //int numberIntegerInfeasibilities=0; // without odd ones
953    numberStrongDone=0;
954    numberUnfinished=0;
955    numberStrongInfeasible=0;
956    numberStrongIterations=0;
957   
958    // We may go round this loop twice (only if we think we have solution)
959    for (int iPass=0;iPass<2;iPass++) {
960     
961      // compute current state
962      //int numberObjectInfeasibilities; // just odd ones
963      //model->feasibleSolution(
964      //                      numberIntegerInfeasibilities,
965      //                      numberObjectInfeasibilities);
966      const double * hotstartSolution = model->hotstartSolution();
967      const int * hotstartPriorities = model->hotstartPriorities();
968     
969      // Some objects may compute an estimate of best solution from here
970      estimatedDegradation=0.0; 
971      numberUnsatisfied_ = 0;
972      int bestPriority=INT_MAX;
973      /*
974        Scan for branching objects that indicate infeasibility. Choose the best
975        maximumStrong candidates, using priority as the first criteria, then
976        integer infeasibility.
977       
978        The algorithm is to fill the choice array with a set of good candidates (by
979        infeasibility) with priority bestPriority.  Finding a candidate with
980        priority better (less) than bestPriority flushes the choice array. (This
981        serves as initialization when the first candidate is found.)
982       
983        A new candidate is added to choices only if its infeasibility exceeds the
984        current max infeasibility (mostAway). When a candidate is added, it
985        replaces the candidate with the smallest infeasibility (tracked by
986        iSmallest).
987      */
988      int iSmallest = 0;
989      double mostAway = 1.0e-100;
990      for (i = 0 ; i < maximumStrong ; i++)
991        choice[i].possibleBranch = NULL ;
992      numberStrong=0;
993      bool canDoOneHot=false;
994      for (i=0;i<numberObjects;i++) {
995        CbcObject * object = model->modifiableObject(i);
996        int preferredWay;
997        double infeasibility = object->infeasibility(preferredWay);
998        int priorityLevel = object->priority();
999        if (hotstartSolution) {
1000          // we are doing hot start
1001          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1002          if (thisOne) {
1003            int iColumn = thisOne->modelSequence();
1004            bool canDoThisHot=true;
1005            double targetValue = hotstartSolution[iColumn];
1006            if (saveUpper[iColumn]>saveLower[iColumn]) {
1007              double value = saveSolution[iColumn];
1008              if (hotstartPriorities)
1009                priorityLevel=hotstartPriorities[iColumn]; 
1010              //double originalLower = thisOne->originalLower();
1011              //double originalUpper = thisOne->originalUpper();
1012              // switch off if not possible
1013              if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
1014                /* priority outranks rest always if negative
1015                   otherwise can be downgraded if at correct level.
1016                   Infeasibility may be increased to choose 1.0 values first.
1017                   choose one near wanted value
1018                */
1019                if (fabs(value-targetValue)>integerTolerance) {
1020                  infeasibility = 1.0-fabs(value-targetValue);
1021                  if (targetValue==1.0)
1022                    infeasibility += 1.0;
1023                  if (value>targetValue) {
1024                    preferredWay=-1;
1025                  } else {
1026                    preferredWay=1;
1027                  }
1028                  priorityLevel = CoinAbs(priorityLevel);
1029                } else if (priorityLevel<0) {
1030                  priorityLevel = CoinAbs(priorityLevel);
1031                  if (targetValue==saveLower[iColumn]) {
1032                    infeasibility = integerTolerance+1.0e-12;
1033                    preferredWay=-1;
1034                  } else if (targetValue==saveUpper[iColumn]) {
1035                    infeasibility = integerTolerance+1.0e-12;
1036                    preferredWay=1;
1037                  } else {
1038                    // can't
1039                    priorityLevel += 10000000;
1040                    canDoThisHot=false;
1041                  }
1042                } else {
1043                  priorityLevel += 10000000;
1044                  canDoThisHot=false;
1045                }
1046              } else {
1047                // switch off if not possible
1048                canDoThisHot=false;
1049              }
1050              if (canDoThisHot)
1051                canDoOneHot=true;
1052            } else if (targetValue<saveLower[iColumn]||targetValue>saveUpper[iColumn]) {
1053            }
1054          } else {
1055            priorityLevel += 10000000;
1056          }
1057        }
1058        if (infeasibility) {
1059          // Increase estimated degradation to solution
1060          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
1061          numberUnsatisfied_++;
1062          // Better priority? Flush choices.
1063          if (priorityLevel<bestPriority) {
1064            int j;
1065            iSmallest=0;
1066            for (j=0;j<maximumStrong;j++) {
1067              choice[j].upMovement=0.0;
1068              delete choice[j].possibleBranch;
1069              choice[j].possibleBranch=NULL;
1070            }
1071            bestPriority = priorityLevel;
1072            mostAway=1.0e-100;
1073            numberStrong=0;
1074          } else if (priorityLevel>bestPriority) {
1075            continue;
1076          }
1077          // Check for suitability based on infeasibility.
1078          if (infeasibility>mostAway) {
1079            //add to list
1080            choice[iSmallest].upMovement=infeasibility;
1081            delete choice[iSmallest].possibleBranch;
1082            choice[iSmallest].possibleBranch=object->createBranch(preferredWay);
1083            numberStrong = CoinMax(numberStrong,iSmallest+1);
1084            // Save which object it was
1085            choice[iSmallest].objectNumber=i;
1086            int j;
1087            iSmallest=-1;
1088            mostAway = 1.0e50;
1089            for (j=0;j<maximumStrong;j++) {
1090              if (choice[j].upMovement<mostAway) {
1091                mostAway=choice[j].upMovement;
1092                iSmallest=j;
1093              }
1094            }
1095          }
1096        }
1097      }
1098      if (!canDoOneHot&&hotstartSolution) {
1099        // switch off as not possible
1100        hotstartSolution=NULL;
1101        model->setHotstartSolution(NULL,NULL);
1102      }
1103      if (numberUnsatisfied_) {
1104        // some infeasibilities - go to next steps
1105        break;
1106      } else if (!iPass) {
1107        // looks like a solution - get paranoid
1108        bool roundAgain=false;
1109        // get basis
1110        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1111        if (!ws)
1112          break;
1113        for (i=0;i<numberColumns;i++) {
1114          double value = saveSolution[i];
1115          if (value<lower[i]) {
1116            saveSolution[i]=lower[i];
1117            roundAgain=true;
1118            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1119          } else if (value>upper[i]) {
1120            saveSolution[i]=upper[i];
1121            roundAgain=true;
1122            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1123          } 
1124        }
1125        if (roundAgain&&saveNumberStrong) {
1126          // restore basis
1127          solver->setWarmStart(ws);
1128          delete ws;
1129          solver->resolve();
1130          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1131          model->reserveCurrentSolution(saveSolution);
1132          if (!solver->isProvenOptimal()) {
1133            // infeasible
1134            anyAction=-2;
1135            break;
1136          }
1137        } else {
1138          delete ws;
1139          break;
1140        }
1141      }
1142    }
1143    /* Some solvers can do the strong branching calculations faster if
1144       they do them all at once.  At present only Clp does for ordinary
1145       integers but I think this coding would be easy to modify
1146    */
1147    bool allNormal=true; // to say if we can do fast strong branching
1148    // Say which one will be best
1149    int bestChoice=0;
1150    double worstInfeasibility=0.0;
1151    for (i=0;i<numberStrong;i++) {
1152      choice[i].numIntInfeasUp = numberUnsatisfied_;
1153      choice[i].numIntInfeasDown = numberUnsatisfied_;
1154      choice[i].fix=0; // say not fixed
1155      if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1156        allNormal=false; // Something odd so lets skip clever fast branching
1157      if ( !model->object(choice[i].objectNumber)->boundBranch())
1158        numberStrong=0; // switch off
1159      if ( choice[i].possibleBranch->numberBranches()>2)
1160        numberStrong=0; // switch off
1161      // Do best choice in case switched off
1162      if (choice[i].upMovement>worstInfeasibility) {
1163        worstInfeasibility=choice[i].upMovement;
1164        bestChoice=i;
1165      }
1166    }
1167    // If we have hit max time don't do strong branching
1168    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1169                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1170    // also give up if we are looping round too much
1171    if (hitMaxTime||numberPassesLeft<=0)
1172      numberStrong=0;
1173    /*
1174      Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1175      fall through to simple branching.
1176     
1177      Setup for strong branching involves saving the current basis (for restoration
1178      afterwards) and setting up for hot starts.
1179    */
1180    if (numberStrong&&saveNumberStrong) {
1181     
1182      bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1183      solveAll=true;
1184      // worth trying if too many times
1185      // Save basis
1186      CoinWarmStart * ws = solver->getWarmStart();
1187      // save limit
1188      int saveLimit;
1189      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1190      if (beforeSolution&&saveLimit<100)
1191        solver->setIntParam(OsiMaxNumIterationHotStart,100); // go to end
1192#     ifdef COIN_HAS_CLP     
1193      /* If we are doing all strong branching in one go then we create new arrays
1194         to store information.  If clp NULL then doing old way.
1195         Going down -
1196         outputSolution[2*i] is final solution.
1197         outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1198         outputStuff[2*i+numberStrong] is number iterations
1199         On entry newUpper[i] is new upper bound, on exit obj change
1200         Going up -
1201         outputSolution[2*i+1] is final solution.
1202         outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1203         outputStuff[2*i+1+numberStrong] is number iterations
1204       On entry newLower[i] is new lower bound, on exit obj change
1205      */
1206      OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1207      ClpSimplex * clp=NULL;
1208      double * newLower = NULL;
1209      double * newUpper = NULL;
1210      double ** outputSolution=NULL;
1211      int * outputStuff=NULL;
1212      // Go back to normal way if user wants it
1213      if (osiclp&&(osiclp->specialOptions()&16)!=0&&osiclp->specialOptions()>0)
1214        allNormal=false;
1215      if (osiclp&&!allNormal) {
1216        // say do fast
1217        int easy=1;
1218        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1219      }
1220      if (osiclp&& allNormal) {
1221        clp = osiclp->getModelPtr();
1222        // Clp - do a different way
1223        newLower = new double[numberStrong];
1224        newUpper = new double[numberStrong];
1225        outputSolution = new double * [2*numberStrong];
1226        outputStuff = new int [4*numberStrong];
1227        int * which = new int[numberStrong];
1228        int startFinishOptions;
1229        int specialOptions = osiclp->specialOptions();
1230        int clpOptions = clp->specialOptions();
1231        int returnCode=0;
1232#define CRUNCH
1233#ifdef CRUNCH
1234        // Crunch down problem
1235        int numberRows = clp->numberRows();
1236        // Use dual region
1237        double * rhs = clp->dualRowSolution();
1238        int * whichRow = new int[3*numberRows];
1239        int * whichColumn = new int[2*numberColumns];
1240        int nBound;
1241        ClpSimplex * small = ((ClpSimplexOther *) clp)->crunch(rhs,whichRow,whichColumn,nBound,true);
1242        if (!small) {
1243          anyAction=-2;
1244          //printf("XXXX Inf by inspection\n");
1245          delete [] whichColumn;
1246          whichColumn=NULL;
1247          delete [] whichRow;
1248          whichRow=NULL;
1249          break;
1250        } else {
1251          clp = small;
1252        }
1253#else
1254        int saveLogLevel = clp->logLevel();
1255        int saveMaxIts = clp->maximumIterations();
1256#endif
1257        clp->setLogLevel(0);
1258        if((specialOptions&1)==0) {
1259          startFinishOptions=0;
1260          clp->setSpecialOptions(clpOptions|(64|1024));
1261        } else {
1262          startFinishOptions=1+2+4;
1263          //startFinishOptions=1+4; // for moment re-factorize
1264          if((specialOptions&4)==0) 
1265            clp->setSpecialOptions(clpOptions|(64|128|512|1024|4096));
1266          else
1267            clp->setSpecialOptions(clpOptions|(64|128|512|1024|2048|4096));
1268        }
1269        // User may want to clean up before strong branching
1270        if ((clp->specialOptions()&32)!=0) {
1271          clp->primal(1);
1272          if (clp->numberIterations())
1273            model->messageHandler()->message(CBC_ITERATE_STRONG,model->messages())
1274              << clp->numberIterations()
1275              <<CoinMessageEol;
1276        }
1277        clp->setMaximumIterations(saveLimit);
1278#ifdef CRUNCH
1279        int * backColumn = whichColumn+numberColumns;
1280#endif
1281        for (i=0;i<numberStrong;i++) {
1282          int iObject = choice[i].objectNumber;
1283          const CbcObject * object = model->object(iObject);
1284          const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1285          int iSequence = simple->modelSequence();
1286          newLower[i]= ceil(saveSolution[iSequence]);
1287          newUpper[i]= floor(saveSolution[iSequence]);
1288#ifdef CRUNCH
1289          iSequence = backColumn[iSequence];
1290          assert (iSequence>=0);
1291#endif
1292          which[i]=iSequence;
1293          outputSolution[2*i]= new double [numberColumns];
1294          outputSolution[2*i+1]= new double [numberColumns];
1295        }
1296        //clp->writeMps("bad");
1297        returnCode=clp->strongBranching(numberStrong,which,
1298                                            newLower, newUpper,outputSolution,
1299                                            outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1300                                            startFinishOptions);
1301#ifndef CRUNCH
1302        clp->setSpecialOptions(clpOptions); // restore
1303        clp->setMaximumIterations(saveMaxIts);
1304        clp->setLogLevel(saveLogLevel);
1305#endif
1306        if (returnCode==-2) {
1307          // bad factorization!!!
1308          // Doing normal way
1309          // Mark hot start
1310          solver->markHotStart();
1311          clp = NULL;
1312        } else {
1313#ifdef CRUNCH
1314          // extract solution
1315          //bool checkSol=true;
1316          for (i=0;i<numberStrong;i++) {
1317            int iObject = choice[i].objectNumber;
1318            const CbcObject * object = model->object(iObject);
1319            const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1320            int iSequence = simple->modelSequence();
1321            which[i]=iSequence;
1322            double * sol = outputSolution[2*i];
1323            double * sol2 = outputSolution[2*i+1];
1324            //bool x=true;
1325            //bool x2=true;
1326            for (int iColumn=numberColumns-1;iColumn>=0;iColumn--) {
1327              int jColumn = backColumn[iColumn];
1328              if (jColumn>=0) {
1329                sol[iColumn]=sol[jColumn];
1330                sol2[iColumn]=sol2[jColumn];
1331              } else {
1332                sol[iColumn]=saveSolution[iColumn];
1333                sol2[iColumn]=saveSolution[iColumn];
1334              }
1335            }
1336          }
1337#endif
1338        }
1339#ifdef CRUNCH
1340        delete [] whichColumn;
1341        delete [] whichRow;
1342        delete small;
1343#endif
1344        delete [] which;
1345      } else {
1346        // Doing normal way
1347        // Mark hot start
1348        solver->markHotStart();
1349      }
1350#     else      /* COIN_HAS_CLP */
1351
1352      OsiSolverInterface *clp = NULL ;
1353      double **outputSolution = NULL ;
1354      int *outputStuff = NULL ;
1355      double * newLower = NULL ;
1356      double * newUpper = NULL ;
1357
1358      solver->markHotStart();
1359
1360#     endif     /* COIN_HAS_CLP */
1361      /*
1362        Open a loop to do the strong branching LPs. For each candidate variable,
1363        solve an LP with the variable forced down, then up. If a direction turns
1364        out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1365        force the objective change to be big (1.0e100). If we determine the problem
1366        is infeasible, or find a monotone variable, escape the loop.
1367       
1368        TODO: The `restore bounds' part might be better encapsulated as an
1369        unbranch() method. Branching objects more exotic than simple integers
1370        or cliques might not restrict themselves to variable bounds.
1371
1372        TODO: Virtuous solvers invalidate the current solution (or give bogus
1373        results :-) when the bounds are changed out from under them. So we
1374        need to do all the work associated with finding a new solution before
1375        restoring the bounds.
1376      */
1377      for (i = 0 ; i < numberStrong ; i++)
1378        { double objectiveChange ;
1379        double newObjectiveValue=1.0e100;
1380        // status is 0 finished, 1 infeasible and other
1381        int iStatus;
1382        /*
1383          Try the down direction first. (Specify the initial branching alternative as
1384          down with a call to way(-1). Each subsequent call to branch() performs the
1385          specified branch and advances the branch object state to the next branch
1386          alternative.)
1387        */
1388        if (!clp) {
1389          choice[i].possibleBranch->way(-1) ;
1390          choice[i].possibleBranch->branch() ;
1391          bool feasible=true;
1392          if (checkFeasibility) {
1393            // check branching did not make infeasible
1394            int iColumn;
1395            int numberColumns = solver->getNumCols();
1396            const double * columnLower = solver->getColLower();
1397            const double * columnUpper = solver->getColUpper();
1398            for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1399              if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1400                feasible=false;
1401            }
1402          }
1403          if (feasible) {
1404            solver->solveFromHotStart() ;
1405            numberStrongDone++;
1406            numberStrongIterations += solver->getIterationCount();
1407            /*
1408              We now have an estimate of objective degradation that we can use for strong
1409              branching. If we're over the cutoff, the variable is monotone up.
1410              If we actually made it to optimality, check for a solution, and if we have
1411              a good one, call setBestSolution to process it. Note that this may reduce the
1412              cutoff, so we check again to see if we can declare this variable monotone.
1413            */
1414            if (solver->isProvenOptimal())
1415              iStatus=0; // optimal
1416            else if (solver->isIterationLimitReached()
1417                     &&!solver->isDualObjectiveLimitReached())
1418              iStatus=2; // unknown
1419            else
1420              iStatus=1; // infeasible
1421            newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1422            choice[i].numItersDown = solver->getIterationCount();
1423          } else {
1424            iStatus=1; // infeasible
1425            newObjectiveValue = 1.0e100;
1426            choice[i].numItersDown = 0;
1427          }
1428        } else {
1429          iStatus = outputStuff[2*i];
1430          choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1431          numberStrongDone++;
1432          numberStrongIterations += choice[i].numItersDown;
1433          newObjectiveValue = objectiveValue+newUpper[i];
1434          solver->setColSolution(outputSolution[2*i]);
1435        }
1436        objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
1437        if (!iStatus) {
1438          choice[i].finishedDown = true ;
1439          if (newObjectiveValue>=model->getCutoff()) {
1440            objectiveChange = 1.0e100; // say infeasible
1441            numberStrongInfeasible++;
1442          } else {
1443            // See if integer solution
1444            if (model->feasibleSolution(choice[i].numIntInfeasDown,
1445                                        choice[i].numObjInfeasDown)
1446                &&model->problemFeasibility()->feasible(model,-1)>=0) {
1447              model->setBestSolution(CBC_STRONGSOL,
1448                                     newObjectiveValue,
1449                                     solver->getColSolution()) ;
1450              // only needed for odd solvers
1451              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1452              objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
1453              model->setLastHeuristic(NULL);
1454              model->incrementUsed(solver->getColSolution());
1455              if (newObjectiveValue >= model->getCutoff()) {    //  *new* cutoff
1456                objectiveChange = 1.0e100 ;
1457                numberStrongInfeasible++;
1458              }
1459            }
1460          }
1461        } else if (iStatus==1) {
1462          objectiveChange = 1.0e100 ;
1463          numberStrongInfeasible++;
1464        } else {
1465          // Can't say much as we did not finish
1466          choice[i].finishedDown = false ;
1467          numberUnfinished++;
1468        }
1469        choice[i].downMovement = objectiveChange ;
1470       
1471        // restore bounds
1472        if (!clp)
1473          { for (int j=0;j<numberColumns;j++) {
1474            if (saveLower[j] != lower[j])
1475              solver->setColLower(j,saveLower[j]);
1476            if (saveUpper[j] != upper[j])
1477              solver->setColUpper(j,saveUpper[j]);
1478          }
1479          }
1480        //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1481        //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
1482        //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
1483        //     choice[i].numObjInfeasDown);
1484       
1485        // repeat the whole exercise, forcing the variable up
1486        if (!clp) {
1487          bool feasible=true;
1488          // If odd branching then maybe just one possibility
1489          if(choice[i].possibleBranch->numberBranchesLeft()>0) {
1490            choice[i].possibleBranch->branch();
1491            if (checkFeasibility) {
1492              // check branching did not make infeasible
1493              int iColumn;
1494              int numberColumns = solver->getNumCols();
1495              const double * columnLower = solver->getColLower();
1496              const double * columnUpper = solver->getColUpper();
1497              for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1498                if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1499                  feasible=false;
1500              }
1501            }
1502          } else {
1503            // second branch infeasible
1504            feasible=false;
1505          }
1506          if (feasible) {
1507            solver->solveFromHotStart() ;
1508            numberStrongDone++;
1509            numberStrongIterations += solver->getIterationCount();
1510            /*
1511              We now have an estimate of objective degradation that we can use for strong
1512              branching. If we're over the cutoff, the variable is monotone up.
1513              If we actually made it to optimality, check for a solution, and if we have
1514              a good one, call setBestSolution to process it. Note that this may reduce the
1515              cutoff, so we check again to see if we can declare this variable monotone.
1516            */
1517            if (solver->isProvenOptimal())
1518              iStatus=0; // optimal
1519            else if (solver->isIterationLimitReached()
1520                     &&!solver->isDualObjectiveLimitReached())
1521              iStatus=2; // unknown
1522            else
1523              iStatus=1; // infeasible
1524            newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1525            choice[i].numItersUp = solver->getIterationCount();
1526          } else {
1527            iStatus=1; // infeasible
1528            newObjectiveValue = 1.0e100;
1529            choice[i].numItersDown = 0;
1530          }
1531        } else {
1532          iStatus = outputStuff[2*i+1];
1533          choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
1534          numberStrongDone++;
1535          numberStrongIterations += choice[i].numItersUp;
1536          newObjectiveValue = objectiveValue+newLower[i];
1537          solver->setColSolution(outputSolution[2*i+1]);
1538        }
1539        objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
1540        if (!iStatus) {
1541          choice[i].finishedUp = true ;
1542          if (newObjectiveValue>=model->getCutoff()) {
1543            objectiveChange = 1.0e100; // say infeasible
1544            numberStrongInfeasible++;
1545          } else {
1546            // See if integer solution
1547            if (model->feasibleSolution(choice[i].numIntInfeasUp,
1548                                        choice[i].numObjInfeasUp)
1549                &&model->problemFeasibility()->feasible(model,-1)>=0) {
1550              model->setBestSolution(CBC_STRONGSOL,
1551                                     newObjectiveValue,
1552                                     solver->getColSolution()) ;
1553              // only needed for odd solvers
1554              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1555              objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
1556              model->setLastHeuristic(NULL);
1557              model->incrementUsed(solver->getColSolution());
1558              if (newObjectiveValue >= model->getCutoff()) {    //  *new* cutoff
1559                objectiveChange = 1.0e100 ;
1560                numberStrongInfeasible++;
1561              }
1562            }
1563          }
1564        } else if (iStatus==1) {
1565          objectiveChange = 1.0e100 ;
1566          numberStrongInfeasible++;
1567        } else {
1568          // Can't say much as we did not finish
1569          choice[i].finishedUp = false ;
1570          numberUnfinished++;
1571        }
1572        choice[i].upMovement = objectiveChange ;
1573       
1574        // restore bounds
1575        if (!clp)
1576          { for (int j=0;j<numberColumns;j++) {
1577            if (saveLower[j] != lower[j])
1578              solver->setColLower(j,saveLower[j]);
1579            if (saveUpper[j] != upper[j])
1580              solver->setColUpper(j,saveUpper[j]);
1581          }
1582          }
1583       
1584        //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1585        //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
1586        //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
1587        //     choice[i].numObjInfeasUp);
1588       
1589        /*
1590          End of evaluation for this candidate variable. Possibilities are:
1591          * Both sides below cutoff; this variable is a candidate for branching.
1592          * Both sides infeasible or above the objective cutoff: no further action
1593          here. Break from the evaluation loop and assume the node will be purged
1594          by the caller.
1595          * One side below cutoff: Install the branch (i.e., fix the variable). Break
1596          from the evaluation loop and assume the node will be reoptimised by the
1597          caller.
1598        */
1599        // reset
1600        choice[i].possibleBranch->resetNumberBranchesLeft();
1601        if (choice[i].upMovement<1.0e100) {
1602          if(choice[i].downMovement<1.0e100) {
1603            // feasible - no action
1604          } else {
1605            // up feasible, down infeasible
1606            anyAction=-1;
1607            //printf("Down infeasible for choice %d sequence %d\n",i,
1608            // model->object(choice[i].objectNumber)->columnNumber());
1609            if (!solveAll) {
1610              choice[i].possibleBranch->way(1);
1611              choice[i].possibleBranch->branch();
1612              break;
1613            } else {
1614              choice[i].fix=1;
1615            }
1616          }
1617        } else {
1618          if(choice[i].downMovement<1.0e100) {
1619            // down feasible, up infeasible
1620            anyAction=-1;
1621            //printf("Up infeasible for choice %d sequence %d\n",i,
1622            // model->object(choice[i].objectNumber)->columnNumber());
1623            if (!solveAll) {
1624              choice[i].possibleBranch->way(-1);
1625              choice[i].possibleBranch->branch();
1626              break;
1627            } else {
1628              choice[i].fix=-1;
1629            }
1630          } else {
1631            // neither side feasible
1632            anyAction=-2;
1633            //printf("Both infeasible for choice %d sequence %d\n",i,
1634            // model->object(choice[i].objectNumber)->columnNumber());
1635            break;
1636          }
1637        }
1638        bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1639                            model->getDblParam(CbcModel::CbcMaximumSeconds));
1640        if (hitMaxTime) {
1641          numberStrong=i+1;
1642          break;
1643        }
1644        }
1645      if (!clp) {
1646        // Delete the snapshot
1647        solver->unmarkHotStart();
1648      } else {
1649        delete [] newLower;
1650        delete [] newUpper;
1651        delete [] outputStuff;
1652        int i;
1653        for (i=0;i<2*numberStrong;i++)
1654          delete [] outputSolution[i];
1655        delete [] outputSolution;
1656      }
1657      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
1658      // restore basis
1659      solver->setWarmStart(ws);
1660      // Unless infeasible we will carry on
1661      // But we could fix anyway
1662      if (anyAction==-1&&solveAll) {
1663        // apply and take off
1664        for (i = 0 ; i < numberStrong ; i++) {
1665          if (choice[i].fix) {
1666            choice[i].possibleBranch->way(choice[i].fix) ;
1667            choice[i].possibleBranch->branch() ;
1668          }
1669        }
1670        bool feasible=true;
1671        if (checkFeasibility) {
1672          // check branching did not make infeasible
1673          int iColumn;
1674          int numberColumns = solver->getNumCols();
1675          const double * columnLower = solver->getColLower();
1676          const double * columnUpper = solver->getColUpper();
1677          for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1678            if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1679              feasible=false;
1680          }
1681        }
1682        if (feasible) {
1683          // can do quick optimality check
1684          int easy=2;
1685          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1686          solver->resolve() ;
1687          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1688          feasible = solver->isProvenOptimal();
1689        }
1690        if (feasible) {
1691          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1692          model->reserveCurrentSolution(saveSolution);
1693          memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
1694          memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
1695          // Clean up all candidates whih are fixed
1696          int numberLeft=0;
1697          for (i = 0 ; i < numberStrong ; i++) {
1698            CbcStrongInfo thisChoice = choice[i];
1699            choice[i].possibleBranch=NULL;
1700            const CbcObject * object = model->object(thisChoice.objectNumber);
1701            int preferredWay;
1702            double infeasibility = object->infeasibility(preferredWay);
1703            if (!infeasibility) {
1704              // take out
1705              delete thisChoice.possibleBranch;
1706            } else {
1707              choice[numberLeft++]=thisChoice;
1708            }
1709          }
1710          numberStrong=numberLeft;
1711          for (;i<maximumStrong;i++) {
1712            delete choice[i].possibleBranch;
1713            choice[i].possibleBranch=NULL;
1714          }
1715          // If all fixed then round again
1716          if (!numberLeft) {
1717            finished=false;
1718            numberStrong=0;
1719            saveNumberStrong=0;
1720            maximumStrong=1;
1721          } else {
1722            anyAction=0;
1723          }
1724          // If these two uncommented then different action
1725          anyAction=-1;
1726          finished=true;
1727          //printf("some fixed but continuing %d left\n",numberLeft);
1728        } else {
1729          anyAction=-2; // say infeasible
1730        }
1731      }
1732      delete ws;
1733      int numberNodes = model->getNodeCount();
1734      // update number of strong iterations etc
1735      model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
1736                                 anyAction==-2 ? 0:numberStrongInfeasible,anyAction==-2);
1737     
1738      /*
1739        anyAction >= 0 indicates that strong branching didn't produce any monotone
1740        variables. Sift through the candidates for the best one.
1741       
1742        QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
1743        local to this code block. Perhaps should be numberNodes_ from model?
1744        Unclear what this calculation is doing.
1745      */
1746      if (anyAction>=0) {
1747       
1748        // get average cost per iteration and assume stopped ones
1749        // would stop after 50% more iterations at average cost??? !!! ???
1750        double averageCostPerIteration=0.0;
1751        double totalNumberIterations=1.0;
1752        int smallestNumberInfeasibilities=INT_MAX;
1753        for (i=0;i<numberStrong;i++) {
1754          totalNumberIterations += choice[i].numItersDown +
1755            choice[i].numItersUp ;
1756          averageCostPerIteration += choice[i].downMovement +
1757            choice[i].upMovement;
1758          smallestNumberInfeasibilities= 
1759            CoinMin(CoinMin(choice[i].numIntInfeasDown ,
1760                            choice[i].numIntInfeasUp ),
1761                    smallestNumberInfeasibilities);
1762        }
1763        //if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
1764        //numberNodes=1000000; // switch off search for better solution
1765        numberNodes=1000000; // switch off anyway
1766        averageCostPerIteration /= totalNumberIterations;
1767        // all feasible - choose best bet
1768       
1769        // New method does all at once so it can be more sophisticated
1770        // in deciding how to balance actions.
1771        // But it does need arrays
1772        double * changeUp = new double [numberStrong];
1773        int * numberInfeasibilitiesUp = new int [numberStrong];
1774        double * changeDown = new double [numberStrong];
1775        int * numberInfeasibilitiesDown = new int [numberStrong];
1776        CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
1777        for (i = 0 ; i < numberStrong ; i++) {
1778          int iColumn = choice[i].possibleBranch->variable() ;
1779          model->messageHandler()->message(CBC_STRONG,model->messages())
1780            << i << iColumn
1781            <<choice[i].downMovement<<choice[i].numIntInfeasDown
1782            <<choice[i].upMovement<<choice[i].numIntInfeasUp
1783            <<choice[i].possibleBranch->value()
1784            <<CoinMessageEol;
1785          changeUp[i]=choice[i].upMovement;
1786          numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1787          changeDown[i]=choice[i].downMovement;
1788          numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1789          objects[i] = choice[i].possibleBranch;
1790        }
1791        int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
1792                                               changeUp,numberInfeasibilitiesUp,
1793                                               changeDown,numberInfeasibilitiesDown,
1794                                               objectiveValue_);
1795        // move branching object and make sure it will not be deleted
1796        if (whichObject>=0) {
1797          branch_ = objects[whichObject];
1798          if (model->messageHandler()->logLevel()>3) 
1799            printf("Choosing column %d\n",choice[whichObject].possibleBranch->variable()) ;
1800          choice[whichObject].possibleBranch=NULL;
1801        }
1802        delete [] changeUp;
1803        delete [] numberInfeasibilitiesUp;
1804        delete [] changeDown;
1805        delete [] numberInfeasibilitiesDown;
1806        delete [] objects;
1807      }
1808#     ifdef COIN_HAS_CLP
1809      if (osiclp&&!allNormal) {
1810        // back to normal
1811        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1812      }
1813#     endif
1814    }
1815    /*
1816      Simple branching. Probably just one, but we may have got here
1817      because of an odd branch e.g. a cut
1818    */
1819    else {
1820      // not strong
1821      // C) create branching object
1822      branch_ = choice[bestChoice].possibleBranch;
1823      choice[bestChoice].possibleBranch=NULL;
1824    }
1825  }
1826  // Set guessed solution value
1827  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
1828/*
1829  Cleanup, then we're outta here.
1830*/
1831  if (!model->branchingMethod()||dynamicBranchingObject)
1832    delete decision;
1833   
1834  for (i=0;i<maximumStrong;i++)
1835    delete choice[i].possibleBranch;
1836  delete [] choice;
1837  delete [] saveLower;
1838  delete [] saveUpper;
1839 
1840  // restore solution
1841  solver->setColSolution(saveSolution);
1842  delete [] saveSolution;
1843  return anyAction;
1844}
1845
1846/*
1847  Version for dynamic pseudo costs.
1848
1849  **** For now just return if anything odd
1850  later allow even if odd
1851
1852  The routine scans through the object list of the model looking for objects
1853  that indicate infeasibility. It tests each object using strong branching
1854  and selects the one with the least objective degradation.  A corresponding
1855  branching object is left attached to lastNode.
1856  This version gives preference in evaluation to variables which
1857  have not been evaluated many times.  It also uses numberStrong
1858  to say give up if last few tries have not changed incumbent.
1859  See Achterberg, Koch and Martin.
1860
1861  If strong branching is disabled, a candidate object is chosen essentially
1862  at random (whatever object ends up in pos'n 0 of the candidate array).
1863
1864  If a branching candidate is found to be monotone, bounds are set to fix the
1865  variable and the routine immediately returns (the caller is expected to
1866  reoptimize).
1867
1868  If a branching candidate is found to result in infeasibility in both
1869  directions, the routine immediately returns an indication of infeasibility.
1870
1871  Returns:  0   both branch directions are feasible
1872           -1   branching variable is monotone
1873           -2   infeasible
1874           -3   Use another method
1875
1876           For now just fix on objective from strong branching.
1877*/
1878
1879int CbcNode::chooseDynamicBranch (CbcModel *model, CbcNode *lastNode,
1880                                  OsiSolverBranch * & branches,int numberPassesLeft)
1881
1882{ if (lastNode)
1883    depth_ = lastNode->depth_+1;
1884  else
1885    depth_ = 0;
1886  delete branch_;
1887  branch_=NULL;
1888  OsiSolverInterface * solver = model->solver();
1889  // get information on solver type
1890  const OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
1891  const OsiBabSolver * auxiliaryInfo = dynamic_cast<const OsiBabSolver *> (auxInfo);
1892  if (!auxiliaryInfo) {
1893    // use one from CbcModel
1894    auxiliaryInfo = model->solverCharacteristics();
1895  }
1896  assert (auxiliaryInfo);
1897  //assert(objectiveValue_ == solver->getObjSense()*solver->getObjValue());
1898  double cutoff =model->getCutoff();
1899  double distanceToCutoff=cutoff-objectiveValue_;
1900  const double * lower = solver->getColLower();
1901  const double * upper = solver->getColUpper();
1902  // See what user thinks
1903  int anyAction=model->problemFeasibility()->feasible(model,0);
1904  if (anyAction) {
1905    // will return -2 if infeasible , 0 if treat as integer
1906    return anyAction-1;
1907  }
1908  int i;
1909  int saveStateOfSearch = model->stateOfSearch();
1910  int numberStrong=model->numberStrong();
1911  if (!auxiliaryInfo->warmStart())
1912    numberStrong=0;
1913  // But make more likely to get out after some times
1914  int changeStrategy=numberStrong;
1915  double changeFactor=1.0;
1916  // Use minimum of this and one stored in objects
1917  //int numberBeforeTrust = model->numberBeforeTrust();
1918  int numberObjects = model->numberObjects();
1919  bool checkFeasibility = numberObjects>model->numberIntegers();
1920  // For now return if not simple
1921  if (checkFeasibility)
1922    return -3;
1923  // Return if doing hot start (in BAB sense)
1924  if (model->hotstartSolution()) 
1925    return -3;
1926#define RANGING
1927#ifdef RANGING
1928  // Pass number
1929  int kPass=0;
1930  int numberRows = solver->getNumRows();
1931#endif
1932  int numberColumns = model->getNumCols();
1933  double * saveUpper = new double[numberColumns];
1934  double * saveLower = new double[numberColumns];
1935
1936  // Save solution in case heuristics need good solution later
1937 
1938  double * saveSolution = new double[numberColumns];
1939  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1940  model->reserveCurrentSolution(saveSolution);
1941  /*
1942    Get a branching decision object. Use the default dynamic decision criteria unless
1943    the user has loaded a decision method into the model.
1944  */
1945  CbcBranchDecision *decision = model->branchingMethod();
1946  if (!decision)
1947    decision = new CbcBranchDynamicDecision();
1948  int numberMini=0;
1949  int xPen=0;
1950  int xMark=0;
1951  for (i=0;i<numberColumns;i++) {
1952    saveLower[i] = lower[i];
1953    saveUpper[i] = upper[i];
1954  }
1955  // Get arrays to sort
1956  double * sort = new double[numberObjects];
1957  int * whichObject = new int[numberObjects];
1958  int * objectMark = new int[2*numberObjects+1];
1959  // Arrays with movements
1960  double * upEstimate = new double[numberObjects];
1961  double * downEstimate = new double[numberObjects];
1962  CbcStrongInfo * fixObject = new CbcStrongInfo[numberObjects];
1963  double estimatedDegradation=0.0; 
1964  int numberNodes=model->getNodeCount();
1965  int saveLogLevel = model->logLevel();
1966  if ((numberNodes%500)==0&&false) {
1967    model->setLogLevel(6);
1968    // Get average up and down costs
1969    double averageUp=0.0;
1970    double averageDown=0.0;
1971    int numberUp=0;
1972    int numberDown=0;
1973    int i;
1974    for ( i=0;i<numberObjects;i++) {
1975      CbcObject * object = model->modifiableObject(i);
1976      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
1977        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
1978      assert(dynamicObject);
1979      int  numberUp2=0;
1980      int numberDown2=0;
1981      double up=0.0;
1982      double down=0.0;
1983      if (dynamicObject->numberTimesUp()) {
1984        numberUp++;
1985        averageUp += dynamicObject->upDynamicPseudoCost();
1986        numberUp2 += dynamicObject->numberTimesUp();
1987        up = dynamicObject->upDynamicPseudoCost();
1988      }
1989      if (dynamicObject->numberTimesDown()) {
1990        numberDown++;
1991        averageDown += dynamicObject->downDynamicPseudoCost();
1992        numberDown2 += dynamicObject->numberTimesDown();
1993        down = dynamicObject->downDynamicPseudoCost();
1994      }
1995      if (numberUp2||numberDown2)
1996        printf("col %d - up %d times cost %g, - down %d times cost %g\n",
1997               dynamicObject->columnNumber(),numberUp2,up,numberDown2,down);
1998    }
1999    if (numberUp) 
2000      averageUp /= (double) numberUp;
2001    else
2002      averageUp=1.0;
2003    if (numberDown) 
2004      averageDown /= (double) numberDown;
2005    else
2006      averageDown=1.0;
2007    printf("total - up %d vars average %g, - down %d vars average %g\n",
2008           numberUp,averageUp,numberDown,averageDown);
2009  }
2010  int numberBeforeTrust = model->numberBeforeTrust();
2011  int numberPenalties = model->numberPenalties();
2012  if (numberBeforeTrust>=1000000) {
2013    numberBeforeTrust = numberBeforeTrust % 1000000;
2014    numberPenalties=0;
2015  } else if (numberBeforeTrust<0) {
2016    numberPenalties=numberColumns;
2017    numberBeforeTrust=0;
2018  }
2019  // May go round twice if strong branching fixes all local candidates
2020  bool finished=false;
2021  int numberToFix=0;
2022# ifdef COIN_HAS_CLP
2023  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
2024  int saveClpOptions=0;
2025  if (osiclp) {
2026    // for faster hot start
2027    saveClpOptions = osiclp->specialOptions();
2028    osiclp->setSpecialOptions(saveClpOptions|8192);
2029  }
2030# else
2031  OsiSolverInterface *osiclp = 0 ;
2032# endif
2033  const CglTreeProbingInfo * probingInfo = model->probingInfo();
2034  int saveSearchStrategy2 = model->searchStrategy();
2035  if (saveSearchStrategy2<999) {
2036    // Get average up and down costs
2037    double averageUp=0.0;
2038    double averageDown=0.0;
2039    {
2040      int numberUp=0;
2041      int numberDown=0;
2042      int i;
2043      for ( i=0;i<numberObjects;i++) {
2044        CbcObject * object = model->modifiableObject(i);
2045        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2046          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2047        assert(dynamicObject);
2048        if (dynamicObject->numberTimesUp()) {
2049          numberUp++;
2050          averageUp += dynamicObject->upDynamicPseudoCost();
2051        }
2052        if (dynamicObject->numberTimesDown()) {
2053          numberDown++;
2054          averageDown += dynamicObject->downDynamicPseudoCost();
2055        }
2056      }
2057      if (numberUp) 
2058        averageUp /= (double) numberUp;
2059      else
2060        averageUp=1.0;
2061      if (numberDown) 
2062        averageDown /= (double) numberDown;
2063      else
2064        averageDown=1.0;
2065      for ( i=0;i<numberObjects;i++) {
2066        CbcObject * object = model->modifiableObject(i);
2067        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2068          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2069        assert(dynamicObject);
2070        if (!dynamicObject->numberTimesUp()) 
2071          dynamicObject->setUpDynamicPseudoCost(averageUp);
2072      if (!dynamicObject->numberTimesDown()) 
2073        dynamicObject->setDownDynamicPseudoCost(averageDown);
2074      }
2075    }
2076  } else if (saveSearchStrategy2<1999) {
2077    // pseudo shadow prices
2078    model->pseudoShadow(NULL,NULL);
2079  } else if (saveSearchStrategy2<2999) {
2080    // leave old ones
2081  } else if (saveSearchStrategy2<3999) {
2082    // pseudo shadow prices at root
2083    if (!numberNodes)
2084      model->pseudoShadow(NULL,NULL);
2085  } else {
2086    abort();
2087  }
2088  if (saveSearchStrategy2>=0)
2089    saveSearchStrategy2 = saveSearchStrategy2 % 1000;
2090  if (saveSearchStrategy2==999)
2091    saveSearchStrategy2=-1;
2092  int px[4]={-1,-1,-1,-1};
2093  int saveSearchStrategy = saveSearchStrategy2<99 ? saveSearchStrategy2 : saveSearchStrategy2-100;
2094  bool newWay = saveSearchStrategy2>98;
2095  int numberNotTrusted=0;
2096  int numberStrongDone=0;
2097  int numberUnfinished=0;
2098  int numberStrongInfeasible=0;
2099  int numberStrongIterations=0;
2100  while(!finished) {
2101    finished=true;
2102    decision->initialize(model);
2103    // Some objects may compute an estimate of best solution from here
2104    estimatedDegradation=0.0; 
2105    numberToFix=0;
2106    int numberIntegerInfeasibilities=0; // without odd ones
2107    int numberToDo=0;
2108    int iBestNot=-1;
2109    int iBestGot=-1;
2110    double best=0.0;
2111    numberNotTrusted=0;
2112    numberStrongDone=0;
2113    numberUnfinished=0;
2114    numberStrongInfeasible=0;
2115    numberStrongIterations=0;
2116    int * which = objectMark+numberObjects+1;
2117    int neededPenalties;
2118    int branchingMethod=-1;
2119    // We may go round this loop three times (only if we think we have solution)
2120    for (int iPass=0;iPass<3;iPass++) {
2121     
2122      // compute current state
2123      int numberObjectInfeasibilities; // just odd ones
2124      model->feasibleSolution(
2125                              numberIntegerInfeasibilities,
2126                              numberObjectInfeasibilities);
2127     
2128      // Some objects may compute an estimate of best solution from here
2129      estimatedDegradation=0.0; 
2130      numberUnsatisfied_ = 0;
2131      int bestPriority=INT_MAX;
2132      int number01 = 0;
2133      const fixEntry * entry = NULL;
2134      const int * toZero = NULL;
2135      const int * toOne = NULL;
2136      const int * backward = NULL;
2137      int numberUnsatisProbed=0;
2138      int numberUnsatisNotProbed=0; // 0-1
2139      if (probingInfo) {
2140        number01 = probingInfo->numberIntegers();
2141        entry = probingInfo->fixEntries();
2142        toZero = probingInfo->toZero();
2143        toOne = probingInfo->toOne();
2144        backward = probingInfo->backward();
2145        if (!toZero[number01]) {
2146          // no info
2147          probingInfo=NULL;
2148        }
2149      }
2150      /*
2151        Scan for branching objects that indicate infeasibility. Choose candidates
2152        using priority as the first criteria, then integer infeasibility.
2153       
2154        The algorithm is to fill the array with a set of good candidates (by
2155        infeasibility) with priority bestPriority.  Finding a candidate with
2156        priority better (less) than bestPriority flushes the choice array. (This
2157        serves as initialization when the first candidate is found.)
2158       
2159      */
2160      numberToDo=0;
2161      neededPenalties=0;
2162      iBestNot=-1;
2163      double bestNot=0.0;
2164      iBestGot=-1;
2165      best=0.0;
2166#define PRINT_STUFF -1
2167      for (i=0;i<numberObjects;i++) {
2168        CbcObject * object = model->modifiableObject(i);
2169        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2170          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2171        assert(dynamicObject);
2172        int preferredWay;
2173        double infeasibility = object->infeasibility(preferredWay);
2174        int priorityLevel = object->priority();
2175#define ZERO_ONE 0
2176#define ZERO_FAKE 1.0e20;
2177#if ZERO_ONE==1
2178        // branch on 0-1 first (temp)
2179        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2180          priorityLevel--;
2181#endif
2182#if ZERO_ONE==2
2183        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2184          infeasibility *= ZERO_FAKE;
2185#endif
2186        if (infeasibility) {
2187          // check branching method
2188          if (branchingMethod!=dynamicObject->method()) {
2189            if (branchingMethod==-1)
2190              branchingMethod = dynamicObject->method();
2191            else
2192              branchingMethod = 100;
2193          }
2194          int iColumn = dynamicObject->columnNumber();
2195          if (probingInfo) {
2196            int iSeq = backward[iColumn];
2197            if (iSeq>=0) {
2198              if (toZero[iSeq+1]>toZero[iSeq]) {
2199                numberUnsatisProbed++;
2200              } else {
2201                numberUnsatisNotProbed++;
2202              }
2203            }
2204          }
2205          //double gap = saveUpper[iColumn]-saveLower[iColumn];
2206          // Give precedence to ones with gap of 1.0
2207          //assert(gap>0.0);
2208          //infeasibility /= CoinMin(gap,100.0);
2209          if (!depth_&&false) {
2210            // try closest to 0.5
2211            double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2212            infeasibility = fabs(0.5-part);
2213          }
2214          bool gotDown=false;
2215          int numberThisDown = dynamicObject->numberTimesDown();
2216          if (numberThisDown>=numberBeforeTrust)
2217            gotDown=true;
2218          bool gotUp=false;
2219          int numberThisUp = dynamicObject->numberTimesUp();
2220          if (numberThisUp>=numberBeforeTrust)
2221            gotUp=true;
2222          if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
2223            printf("%d down %d %g up %d %g - infeas %g\n",
2224                   i,numberThisDown,object->downEstimate(),numberThisUp,object->upEstimate(),
2225                   infeasibility);
2226          // Increase estimated degradation to solution
2227          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
2228          downEstimate[i]=object->downEstimate();
2229          upEstimate[i]=object->upEstimate();
2230          numberUnsatisfied_++;
2231          // Better priority? Flush choices.
2232          if (priorityLevel<bestPriority) {
2233            numberToDo=0;
2234            bestPriority = priorityLevel;
2235            iBestGot=-1;
2236            best=0.0;
2237            numberNotTrusted=0;
2238          } else if (priorityLevel>bestPriority) {
2239            continue;
2240          }
2241          if (!gotUp||!gotDown)
2242            numberNotTrusted++;
2243          // Check for suitability based on infeasibility.
2244          if ((gotDown&&gotUp)&&numberStrong>0) {
2245            sort[numberToDo]=-infeasibility;
2246            if (infeasibility>best) {
2247              best=infeasibility;
2248              iBestGot=numberToDo;
2249            }
2250          } else {
2251            objectMark[neededPenalties]=numberToDo;
2252            which[neededPenalties++]=dynamicObject->columnNumber();
2253            int iColumn = dynamicObject->columnNumber();
2254            double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2255            sort[numberToDo]=-10.0*infeasibility;
2256            if (!(numberThisUp+numberThisDown))
2257              sort[numberToDo] *= 100.0; // make even more likely
2258            if (1.0-fabs(part-0.5)>bestNot) {
2259              iBestNot=numberToDo;
2260              bestNot = 1.0-fabs(part-0.5);
2261            }
2262          }
2263          whichObject[numberToDo++]=i;
2264        } else {
2265          // for debug
2266          downEstimate[i]=-1.0;
2267          upEstimate[i]=-1.0;
2268        }
2269      }
2270      if (numberUnsatisfied_) {
2271        if (probingInfo&&false)
2272          printf("nunsat %d, %d probed, %d other 0-1\n",numberUnsatisfied_,
2273                 numberUnsatisProbed,numberUnsatisNotProbed);
2274        // some infeasibilities - go to next steps
2275        break;
2276      } else if (!iPass) {
2277        // may just need resolve
2278        solver->resolve();
2279        memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2280        model->reserveCurrentSolution(saveSolution);
2281        if (!solver->isProvenOptimal()) {
2282          // infeasible
2283          anyAction=-2;
2284          break;
2285        }
2286      } else if (iPass==1) {
2287        // looks like a solution - get paranoid
2288        bool roundAgain=false;
2289        // get basis
2290        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
2291        if (!ws)
2292          break;
2293        double tolerance;
2294        solver->getDblParam(OsiPrimalTolerance,tolerance);
2295        for (i=0;i<numberColumns;i++) {
2296          double value = saveSolution[i];
2297          if (value<lower[i]-tolerance) {
2298            saveSolution[i]=lower[i];
2299            roundAgain=true;
2300            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
2301          } else if (value>upper[i]+tolerance) {
2302            saveSolution[i]=upper[i];
2303            roundAgain=true;
2304            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
2305          } 
2306        }
2307        if (roundAgain) {
2308          // restore basis
2309          solver->setWarmStart(ws);
2310          solver->setColSolution(saveSolution);
2311          delete ws;
2312          bool takeHint;
2313          OsiHintStrength strength;
2314          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
2315          solver->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
2316          solver->resolve();
2317          solver->setHintParam(OsiDoDualInResolve,takeHint,strength) ;
2318          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2319          model->reserveCurrentSolution(saveSolution);
2320          if (!solver->isProvenOptimal()) {
2321            // infeasible
2322            anyAction=-2;
2323            break;
2324          }
2325        } else {
2326          delete ws;
2327          break;
2328        }
2329      }
2330    }
2331    if (anyAction==-2) {
2332      break;
2333    }
2334    bool solveAll=false; // set true to say look at all even if some fixed (experiment)
2335    solveAll=true;
2336    // skip if solution
2337    if (!numberUnsatisfied_)
2338      break;
2339    //bool skipAll = (numberBeforeTrust>20&&numberNodes>20000&&numberNotTrusted==0);
2340    bool skipAll = numberNotTrusted==0;
2341    bool doneHotStart=false;
2342    int searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
2343#ifndef CBC_WEAK_STRONG
2344    if (((numberNodes%20)==0&&searchStrategy!=2)||(model->specialOptions()&8)!=0)
2345      skipAll=false;
2346#endif
2347    if (!newWay) {
2348    // 10 up always use %10, 20 up as 10 and allow penalties
2349    // But adjust depending on ratio of iterations
2350    if (searchStrategy>0&&saveSearchStrategy<10) {
2351      if (numberBeforeTrust>=5&&numberBeforeTrust<=10) {
2352        if (searchStrategy!=2) {
2353          if (depth_>5) {
2354            int numberIterations = model->getIterationCount();
2355            int numberStrongIterations = model->numberStrongIterations();
2356            if (numberStrongIterations>numberIterations+10000) {
2357              searchStrategy=2;
2358              //skipAll=true;
2359            } else if (numberStrongIterations*4+1000<numberIterations||depth_<5) {
2360              searchStrategy=3;
2361              skipAll=false;
2362            }
2363          } else {
2364            searchStrategy=3;
2365            skipAll=false;
2366          }
2367        } else {
2368          //skipAll=true;
2369        }
2370      }
2371    }
2372    } else {
2373    // But adjust depending on ratio of iterations
2374    if (saveSearchStrategy<0) {
2375      // unset
2376      if ((numberNodes%20)==0||(model->specialOptions()&8)!=0) {
2377        // Do numberStrong
2378        searchStrategy=3;
2379      } else if (depth_<5) {
2380        // Do numberStrong
2381        searchStrategy=2;
2382      } else {
2383        int numberIterations = model->getIterationCount();
2384        int numberStrongIterations = model->numberStrongIterations();
2385        int numberRows = solver->getNumRows();
2386        if (numberStrongIterations>numberIterations+CoinMin(10000,10*numberRows)) {
2387          // off
2388          searchStrategy=0;
2389        } else if (numberStrongIterations*4+1000<numberIterations) {
2390          // Do numberStrong if not trusted
2391          searchStrategy=2;
2392        } else {
2393          searchStrategy=1;
2394        }
2395      }
2396    }
2397    if (searchStrategy<3&&(!numberNotTrusted||!searchStrategy))
2398      skipAll=true;
2399    else
2400      skipAll=false;
2401    }
2402    // worth trying if too many times
2403    // Save basis
2404    CoinWarmStart * ws = NULL;
2405    // save limit
2406    int saveLimit=0;
2407    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
2408    if (!skipAll) {
2409      ws = solver->getWarmStart();
2410      int limit=100;
2411#if 0
2412      int averageBranchIterations = model->getIterationCount()/(model->getNodeCount()+1);
2413      if (numberNodes)
2414        limit = CoinMin(CoinMax(limit,2*averageBranchIterations),500);
2415      else
2416        limit = 500;
2417#endif
2418      if ((!saveStateOfSearch||searchStrategy>3)&&saveLimit<limit&&saveLimit==100)
2419        solver->setIntParam(OsiMaxNumIterationHotStart,limit); 
2420    }
2421    // Say which one will be best
2422    int whichChoice=0;
2423    int bestChoice;
2424    if (iBestGot>=0)
2425      bestChoice=iBestGot;
2426    else
2427      bestChoice=iBestNot;
2428    assert (bestChoice>=0);
2429    // If we have hit max time don't do strong branching
2430    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2431                        model->getDblParam(CbcModel::CbcMaximumSeconds));
2432    // also give up if we are looping round too much
2433    if (hitMaxTime||numberPassesLeft<=0||(!numberNotTrusted&&false)||branchingMethod==11) {
2434      int iObject = whichObject[bestChoice];
2435      CbcObject * object = model->modifiableObject(iObject);
2436      int preferredWay;
2437      object->infeasibility(preferredWay);
2438      branch_=object->createBranch(preferredWay);
2439      branch_->way(preferredWay);
2440      delete ws;
2441      ws=NULL;
2442      break;
2443    } else {
2444      // say do fast
2445      int easy=1;
2446      if (!skipAll)
2447        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
2448      int iDo;
2449#ifdef RANGING
2450      if ((skipAll&&numberBeforeTrust&&saveSearchStrategy<20)||saveSearchStrategy<10)
2451        numberPenalties=0;
2452      {
2453        // off penalties if too much
2454        double needed = neededPenalties;
2455        needed *= numberRows;
2456        if (needed>1.0e6&&numberNodes&&saveSearchStrategy<20) {
2457          numberPenalties=0;
2458          neededPenalties=0;
2459        }
2460      }
2461#     ifdef COIN_HAS_CLP
2462      if (osiclp&&numberPenalties&&neededPenalties) {
2463        xPen += neededPenalties;
2464        which--;
2465        which[0]=neededPenalties;
2466        osiclp->passInRanges(which);
2467        // Mark hot start and get ranges
2468        if (kPass) {
2469          // until can work out why solution can go funny
2470          int save = osiclp->specialOptions();
2471          osiclp->setSpecialOptions(save|256);
2472          solver->markHotStart();
2473          osiclp->setSpecialOptions(save);
2474        } else {
2475          solver->markHotStart();
2476        }
2477        assert (auxiliaryInfo->warmStart());
2478        doneHotStart=true;
2479        xMark++;
2480        kPass++;
2481        osiclp->passInRanges(NULL);
2482        const double * downCost=osiclp->upRange();
2483        const double * upCost=osiclp->downRange();
2484        //printf("numberTodo %d needed %d numberpenalties %d\n",numberToDo,neededPenalties,numberPenalties);
2485        double invTrust = 1.0/((double) numberBeforeTrust);
2486        for (int i=0;i<neededPenalties;i++) {
2487          int j = objectMark[i];
2488          int iObject = whichObject[j];
2489          CbcObject * object = model->modifiableObject(iObject);
2490          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2491            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2492          int iSequence=dynamicObject->columnNumber();
2493          double value = saveSolution[iSequence];
2494          value -= floor(value);
2495          double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
2496          double downPenalty = CoinMin(downCost[i],1.0e110)*value;
2497          if (!numberBeforeTrust) {
2498            // override
2499            downEstimate[iObject]=downPenalty;
2500            upEstimate[iObject]=upPenalty;
2501          } else {
2502            int numberThisDown = dynamicObject->numberTimesDown();
2503            if (numberThisDown<numberBeforeTrust) {
2504              double fraction = ((double) numberThisDown)*invTrust;
2505              downEstimate[iObject] = fraction*downEstimate[iObject]+(1.0-fraction)*downPenalty;
2506            }
2507            int numberThisUp = dynamicObject->numberTimesUp();
2508            if (numberThisUp<numberBeforeTrust) {
2509              double fraction = ((double) numberThisUp)*invTrust;
2510              upEstimate[iObject] = fraction*upEstimate[iObject]+(1.0-fraction)*upPenalty;
2511            }
2512          }
2513          sort[j] = - CoinMin(downEstimate[iObject],upEstimate[iObject]);
2514#ifdef CBC_WEAK_STRONG
2515          sort[j] -= 1.0e10; // make more likely to be chosen
2516#endif
2517          //if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
2518          if (!numberNodes)
2519            printf("%d pen down ps %g -> %g up ps %g -> %g\n",
2520                   iObject,downCost[i],downPenalty,upCost[i],upPenalty);
2521        }
2522      } else
2523#     endif     /* COIN_HAS_CLP */
2524      {
2525        if (!skipAll) {
2526          // Mark hot start
2527          solver->markHotStart();
2528          doneHotStart=true;
2529          assert (auxiliaryInfo->warmStart());
2530          xMark++;
2531          //if (solver->isProvenPrimalInfeasible())
2532          //printf("**** Hot start says node infeasible\n");
2533        }
2534        // make sure best will be first
2535        if (iBestGot>=0)
2536          sort[iBestGot]=-1.0e120;
2537      }
2538#else           /* RANGING */
2539      if (!skipAll) {
2540        // Mark hot start
2541        doneHotStart=true;
2542        assert (auxiliaryInfo->warmStart());
2543        solver->markHotStart();
2544        xMark++;
2545      }
2546      // make sure best will be first
2547      if (iBestGot>=0)
2548        sort[iBestGot]=-COIN_DBL_MAX;
2549#endif          /* RANGING */
2550      // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
2551#define ACTION 0
2552#if ACTION<2
2553      if (anyAction)
2554        numberToDo=0; // skip as we will be trying again
2555#endif
2556      // Sort
2557      CoinSort_2(sort,sort+numberToDo,whichObject);
2558      // Change in objective opposite infeasible
2559      double worstFeasible=0.0;
2560      // Just first if strong off
2561      if (!numberStrong)
2562        numberToDo=CoinMin(numberToDo,1);
2563      iDo=0;
2564      int saveLimit2;
2565      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit2);
2566      bool doQuickly = false; // numberToDo>2*numberStrong;
2567      if (searchStrategy==2)
2568        doQuickly=true;
2569      //printf("todo %d, strong %d\n",numberToDo,numberStrong);
2570      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
2571      int numberTest2 = 2*numberStrong;
2572      double distanceToCutoff2 = model->getCutoff()-objectiveValue_;
2573      if (!newWay) {
2574      if (searchStrategy==3) {
2575        // Previously decided we need strong
2576        doQuickly=false;
2577        numberTest = numberStrong;
2578        //numberTest2 = 1000000;
2579      }
2580      if (searchStrategy<0||searchStrategy==1)
2581        //numberTest2 = 1000000;
2582#if 0
2583      if (numberBeforeTrust>20&&(numberNodes>20000||(numberNodes>200&&numberNotTrusted==0))) {
2584        if ((numberNodes%20)!=0) {
2585          numberTest=0;
2586          doQuickly=true;
2587        }
2588      }
2589#else
2590      // Try nearly always off
2591      if (searchStrategy<2) {
2592        if ((numberNodes%20)!=0) {
2593          if ((model->specialOptions()&8)==0) {
2594            numberTest=0;
2595            doQuickly=true;
2596          }
2597        } else {
2598          doQuickly=false;
2599          numberTest=2*numberStrong;
2600          skipAll=false;
2601        }
2602      } else if (searchStrategy!=3) {
2603        doQuickly=true;
2604        numberTest=numberStrong;
2605      }
2606#endif
2607      if (depth_<10&&numberStrong) {
2608        if (searchStrategy!=2) {
2609          doQuickly=false;
2610          int numberRows = solver->getNumRows();
2611          // whether to do this or not is important - think
2612          if (numberRows<300||numberRows+numberColumns<2500) {
2613            if (depth_<7)
2614              numberStrong = CoinMin(3*numberStrong,numberToDo);
2615            if (!depth_) 
2616              numberStrong=CoinMin(6*numberStrong,numberToDo);
2617          }
2618          numberTest=numberStrong;
2619          skipAll=false;
2620        }
2621        model->setStateOfSearch(2); // use min min
2622      }
2623      // could adjust using average iterations per branch
2624      // double average = ((double)model->getIterationCount())/
2625      //((double) model->getNodeCount()+1.0);
2626      // if too many and big then just do 10 its
2627      if (!skipAll&&saveStateOfSearch) {
2628        //if (numberNotTrusted>3*numberStrong&&numberRows>250&&numberColumns>1000&&saveLimit==100)
2629          // off solver->setIntParam(OsiMaxNumIterationHotStart,10);
2630      }
2631      // make negative for test
2632      distanceToCutoff = - distanceToCutoff;
2633      if (numberObjects>-100) {
2634        // larger
2635        distanceToCutoff *= 100.0;
2636      }
2637        distanceToCutoff = -COIN_DBL_MAX;
2638      // Do at least 5 strong
2639      if (numberColumns<1000&&(depth_<15||numberNodes<1000000))
2640        numberTest = CoinMax(numberTest,5);
2641      if ((model->specialOptions()&8)==0) {
2642        if (skipAll) {
2643          numberTest=0;
2644          doQuickly=true;
2645        }
2646      } else {
2647        // do 5 as strong is fixing
2648        numberTest = CoinMax(numberTest,5);
2649      }
2650      } else {
2651      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
2652      int numberTest2 = 2*numberStrong;
2653      if (searchStrategy>=3) {
2654        // Previously decided we need strong
2655        doQuickly=false;
2656        if (depth_<7)
2657          numberStrong *=3;
2658        if (!depth_) 
2659          numberStrong=CoinMin(6*numberStrong,numberToDo);
2660        numberTest = numberStrong;
2661        numberTest2 *= 2;
2662      } else if (searchStrategy==2||(searchStrategy==1&&depth_<6)) {
2663        numberStrong *=2;
2664        if (!depth_) 
2665          numberStrong=CoinMin(2*numberStrong,numberToDo);
2666        numberTest = numberStrong;
2667      } else if (searchStrategy==1&&numberNotTrusted) {
2668        numberTest = numberStrong;
2669      } else {
2670        numberTest=0;
2671        skipAll=true;
2672      }
2673      distanceToCutoff=model->getCutoff()-objectiveValue_;
2674      // make negative for test
2675      distanceToCutoff = - distanceToCutoff;
2676      if (numberObjects>-100) {
2677        // larger
2678        distanceToCutoff *= 100.0;
2679      }
2680      distanceToCutoff = -COIN_DBL_MAX;
2681      if (skipAll) {
2682        numberTest=0;
2683        doQuickly=true;
2684      }
2685      }
2686#if 0
2687      // temp - always switch off
2688      if (0) {
2689        int numberIterations = model->getIterationCount();
2690        int numberStrongIterations = model->numberStrongIterations();
2691        if (numberStrongIterations>numberIterations+10000&&depth_>=5) {
2692          skipAll=true;
2693          newWay=false;
2694          numberTest=0;
2695          doQuickly=true;
2696        }
2697      }
2698      // temp - always switch on
2699      if (0) {
2700        int numberIterations = model->getIterationCount();
2701        int numberStrongIterations = model->numberStrongIterations();
2702        if (2*numberStrongIterations<numberIterations||depth_<=5) {
2703          skipAll=false;
2704          newWay=false;
2705          numberTest=CoinMax(numberTest,numberStrong);
2706          doQuickly=false;
2707        }
2708      }
2709#endif
2710      px[0]=numberTest;
2711      px[1]=numberTest2;
2712      px[2]= doQuickly ? 1 : -1;
2713      px[3]=numberStrong;
2714      if (!newWay) {
2715        if (numberColumns>8*solver->getNumRows()&&false) {
2716          printf("skipAll %c doQuickly %c numberTest %d numberTest2 %d numberNot %d\n",
2717                 skipAll ? 'Y' : 'N',doQuickly ? 'Y' : 'N',numberTest,numberTest2,numberNotTrusted);
2718          numberTest = CoinMin(numberTest,model->numberStrong());
2719          numberTest2 = CoinMin(numberTest2,model->numberStrong());
2720          printf("new test,test2 %d %d\n",numberTest,numberTest2);
2721        }
2722      }
2723      //printf("skipAll %c doQuickly %c numberTest %d numberTest2 %d numberNot %d\n",
2724      //     skipAll ? 'Y' : 'N',doQuickly ? 'Y' : 'N',numberTest,numberTest2,numberNotTrusted);
2725      // See if we want mini tree
2726      bool wantMiniTree=false;
2727      if (model->sizeMiniTree()&&depth_>7&&saveStateOfSearch>0)
2728        wantMiniTree=true;
2729      numberMini=0;
2730      //if (skipAll&&numberTest==0&&doQuickly)
2731      //numberToDo = 1; // trust previous stuff
2732      bool couldChooseFirst = false ; //(skipAll&&numberTest==0&&doQuickly);
2733      //skipAll=false;
2734      for ( iDo=0;iDo<numberToDo;iDo++) {
2735        CbcStrongInfo choice;
2736        int iObject = whichObject[iDo];
2737        CbcObject * object = model->modifiableObject(iObject);
2738        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2739          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2740        int iColumn = dynamicObject->columnNumber();
2741        int preferredWay;
2742        object->infeasibility(preferredWay);
2743        choice.possibleBranch=object->createBranch(preferredWay);
2744        // Save which object it was
2745        choice.objectNumber=iObject;
2746        choice.numIntInfeasUp = numberUnsatisfied_;
2747        choice.numIntInfeasDown = numberUnsatisfied_;
2748        choice.upMovement = upEstimate[iObject];
2749        choice.downMovement = downEstimate[iObject];
2750        assert (choice.upMovement>=0.0);
2751        assert (choice.downMovement>=0.0);
2752        choice.fix=0; // say not fixed
2753        double maxChange = 0.5*(choice.upMovement+choice.downMovement);
2754        maxChange = CoinMin(choice.upMovement,choice.downMovement);
2755        maxChange = CoinMax(choice.upMovement,choice.downMovement);
2756        if (searchStrategy==2)
2757          maxChange = COIN_DBL_MAX;
2758        //maxChange *= 5.0;
2759        if (dynamicObject->method()==1)
2760          maxChange *= 0.1; // probing
2761        // see if can skip strong branching
2762        int canSkip = choice.possibleBranch->fillStrongInfo(choice);
2763        if (!newWay) {
2764          if ((maxChange>distanceToCutoff2)&&(!doQuickly||(numberTest>0&&searchStrategy!=2)))
2765          canSkip=0;
2766        } else {
2767        if (skipAll)
2768          canSkip=1;
2769        else if (numberTest>0&&searchStrategy>=3)
2770          canSkip=0;
2771        }
2772        if (!numberBeforeTrust) {
2773          canSkip=1;
2774        }
2775        if (sort[iDo]<distanceToCutoff)
2776          canSkip=0;
2777        if (((numberTest2<=0&&numberTest<=0)||skipAll)&&sort[iDo]>distanceToCutoff) {
2778          canSkip=1; // always skip
2779          if (iDo>20) {
2780            delete choice.possibleBranch;
2781            choice.possibleBranch=NULL;
2782            break; // give up anyway
2783          }
2784        }
2785        if (model->messageHandler()->logLevel()>3&&numberBeforeTrust) 
2786          dynamicObject->print(1,choice.possibleBranch->value());
2787        // was if (!canSkip)
2788        if (newWay)
2789        numberTest2--;
2790        if (!canSkip) {
2791          //#ifndef RANGING
2792          if (!doneHotStart) {
2793            // Mark hot start
2794            doneHotStart=true;
2795            assert (auxiliaryInfo->warmStart());
2796            solver->markHotStart();
2797            xMark++;
2798          }
2799          //#endif
2800          assert (!couldChooseFirst);
2801          numberTest--;
2802          if (!newWay)
2803          numberTest2--;
2804          // just do a few
2805          //if (canSkip)
2806          //solver->setIntParam(OsiMaxNumIterationHotStart,10);
2807          double objectiveChange ;
2808          double newObjectiveValue=1.0e100;
2809          int j;
2810          // status is 0 finished, 1 infeasible and other
2811          int iStatus;
2812          /*
2813            Try the down direction first. (Specify the initial branching alternative as
2814            down with a call to way(-1). Each subsequent call to branch() performs the
2815            specified branch and advances the branch object state to the next branch
2816            alternative.)
2817          */
2818          choice.possibleBranch->way(-1) ;
2819          decision->saveBranchingObject( choice.possibleBranch);
2820          choice.possibleBranch->branch() ;
2821          solver->solveFromHotStart() ;
2822          bool needHotStartUpdate=false;
2823          numberStrongDone++;
2824          numberStrongIterations += solver->getIterationCount();
2825          /*
2826            We now have an estimate of objective degradation that we can use for strong
2827            branching. If we're over the cutoff, the variable is monotone up.
2828            If we actually made it to optimality, check for a solution, and if we have
2829            a good one, call setBestSolution to process it. Note that this may reduce the
2830            cutoff, so we check again to see if we can declare this variable monotone.
2831          */
2832          if (solver->isProvenOptimal())
2833            iStatus=0; // optimal
2834          else if (solver->isIterationLimitReached()
2835                   &&!solver->isDualObjectiveLimitReached())
2836            iStatus=2; // unknown
2837          else
2838            iStatus=1; // infeasible
2839          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2840          choice.numItersDown = solver->getIterationCount();
2841          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2842          decision->updateInformation( solver,this);
2843          if (!iStatus) {
2844            choice.finishedDown = true ;
2845            if (newObjectiveValue>=cutoff) {
2846              objectiveChange = 1.0e100; // say infeasible
2847              numberStrongInfeasible++;
2848            } else {
2849              // See if integer solution
2850              if (model->feasibleSolution(choice.numIntInfeasDown,
2851                                          choice.numObjInfeasDown)
2852                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2853                if (auxiliaryInfo->solutionAddsCuts()) {
2854                  needHotStartUpdate=true;
2855                  solver->unmarkHotStart();
2856                }
2857                model->setBestSolution(CBC_STRONGSOL,
2858                                       newObjectiveValue,
2859                                       solver->getColSolution()) ;
2860                if (needHotStartUpdate) {
2861                  solver->resolve();
2862                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2863                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2864                  model->feasibleSolution(choice.numIntInfeasDown,
2865                                          choice.numObjInfeasDown);
2866                }
2867                model->setLastHeuristic(NULL);
2868                model->incrementUsed(solver->getColSolution());
2869                cutoff =model->getCutoff();
2870                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
2871                  objectiveChange = 1.0e100 ;
2872                  numberStrongInfeasible++;
2873                }
2874              }
2875            }
2876          } else if (iStatus==1) {
2877            objectiveChange = 1.0e100 ;
2878            numberStrongInfeasible++;
2879          } else {
2880            // Can't say much as we did not finish
2881            choice.finishedDown = false ;
2882            numberUnfinished++;
2883          }
2884          choice.downMovement = objectiveChange ;
2885         
2886          // restore bounds
2887          for ( j=0;j<numberColumns;j++) {
2888            if (saveLower[j] != lower[j])
2889              solver->setColLower(j,saveLower[j]);
2890            if (saveUpper[j] != upper[j])
2891              solver->setColUpper(j,saveUpper[j]);
2892          }
2893          if(needHotStartUpdate) {
2894            needHotStartUpdate = false;
2895            solver->resolve();
2896            //we may again have an integer feasible solution
2897            int numberIntegerInfeasibilities;
2898            int numberObjectInfeasibilities;
2899            if (model->feasibleSolution(
2900                                        numberIntegerInfeasibilities,
2901                                        numberObjectInfeasibilities)) {
2902#ifdef BONMIN
2903              //In this case node has become integer feasible, let us exit the loop
2904              std::cout<<"Node has become integer feasible"<<std::endl;
2905              numberUnsatisfied_ = 0;
2906              break;
2907#endif
2908              double objValue = solver->getObjValue();
2909              model->setBestSolution(CBC_STRONGSOL,
2910                                     objValue,
2911                                     solver->getColSolution()) ;
2912              solver->resolve();
2913              cutoff =model->getCutoff();
2914            }
2915            solver->markHotStart();
2916          }
2917          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2918          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2919          //     choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
2920          //     choice.downMovement,choice.finishedDown,choice.numIntInfeasDown,
2921          //     choice.numObjInfeasDown);
2922         
2923          // repeat the whole exercise, forcing the variable up
2924          decision->saveBranchingObject( choice.possibleBranch);
2925          choice.possibleBranch->branch();
2926          solver->solveFromHotStart() ;
2927          numberStrongDone++;
2928          numberStrongIterations += solver->getIterationCount();
2929          /*
2930            We now have an estimate of objective degradation that we can use for strong
2931            branching. If we're over the cutoff, the variable is monotone up.
2932            If we actually made it to optimality, check for a solution, and if we have
2933            a good one, call setBestSolution to process it. Note that this may reduce the
2934            cutoff, so we check again to see if we can declare this variable monotone.
2935          */
2936          if (solver->isProvenOptimal())
2937            iStatus=0; // optimal
2938          else if (solver->isIterationLimitReached()
2939                   &&!solver->isDualObjectiveLimitReached())
2940            iStatus=2; // unknown
2941          else
2942            iStatus=1; // infeasible
2943          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2944          choice.numItersUp = solver->getIterationCount();
2945          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2946          decision->updateInformation( solver,this);
2947          if (!iStatus) {
2948            choice.finishedUp = true ;
2949            if (newObjectiveValue>=cutoff) {
2950              objectiveChange = 1.0e100; // say infeasible
2951              numberStrongInfeasible++;
2952            } else {
2953              // See if integer solution
2954              if (model->feasibleSolution(choice.numIntInfeasUp,
2955                                          choice.numObjInfeasUp)
2956                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2957#ifdef BONMIN
2958                std::cout<<"Node has become integer feasible"<<std::endl;
2959                numberUnsatisfied_ = 0;
2960                break;
2961#endif
2962                if (auxiliaryInfo->solutionAddsCuts()) {
2963                  needHotStartUpdate=true;
2964                  solver->unmarkHotStart();
2965                }
2966                model->setBestSolution(CBC_STRONGSOL,
2967                                       newObjectiveValue,
2968                                       solver->getColSolution()) ;
2969                if (needHotStartUpdate) {
2970                  solver->resolve();
2971                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2972                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2973                  model->feasibleSolution(choice.numIntInfeasDown,
2974                                          choice.numObjInfeasDown);
2975                }
2976                model->setLastHeuristic(NULL);
2977                model->incrementUsed(solver->getColSolution());
2978                cutoff =model->getCutoff();
2979                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
2980                  objectiveChange = 1.0e100 ;
2981                  numberStrongInfeasible++;
2982                }
2983              }
2984            }
2985          } else if (iStatus==1) {
2986            objectiveChange = 1.0e100 ;
2987            numberStrongInfeasible++;
2988          } else {
2989            // Can't say much as we did not finish
2990            choice.finishedUp = false ;
2991            numberUnfinished++;
2992          }
2993          choice.upMovement = objectiveChange ;
2994         
2995          // restore bounds
2996          for ( j=0;j<numberColumns;j++) {
2997            if (saveLower[j] != lower[j])
2998              solver->setColLower(j,saveLower[j]);
2999            if (saveUpper[j] != upper[j])
3000              solver->setColUpper(j,saveUpper[j]);
3001          }
3002          if(needHotStartUpdate) {
3003            needHotStartUpdate = false;
3004            solver->resolve();
3005            //we may again have an integer feasible solution
3006            int numberIntegerInfeasibilities;
3007            int numberObjectInfeasibilities;
3008            if (model->feasibleSolution(
3009                                        numberIntegerInfeasibilities,
3010                                        numberObjectInfeasibilities)) {
3011              double objValue = solver->getObjValue();
3012              model->setBestSolution(CBC_STRONGSOL,
3013                                     objValue,
3014                                     solver->getColSolution()) ;
3015              solver->resolve();
3016              cutoff =model->getCutoff();
3017            }
3018            solver->markHotStart();
3019          }
3020         
3021          //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3022          //     choice.objectNumber,iStatus,newObjectiveValue,choice.numItersUp,
3023          //     choice.upMovement,choice.finishedUp,choice.numIntInfeasUp,
3024          //     choice.numObjInfeasUp);
3025        }
3026   
3027        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit2); 
3028        /*
3029          End of evaluation for this candidate variable. Possibilities are:
3030          * Both sides below cutoff; this variable is a candidate for branching.
3031          * Both sides infeasible or above the objective cutoff: no further action
3032          here. Break from the evaluation loop and assume the node will be purged
3033          by the caller.
3034          * One side below cutoff: Install the branch (i.e., fix the variable). Break
3035          from the evaluation loop and assume the node will be reoptimised by the
3036          caller.
3037        */
3038        // reset
3039        choice.possibleBranch->resetNumberBranchesLeft();
3040        if (choice.upMovement<1.0e100) {
3041          if(choice.downMovement<1.0e100) {
3042            // In case solution coming in was odd
3043            choice.upMovement = CoinMax(0.0,choice.upMovement);
3044            choice.downMovement = CoinMax(0.0,choice.downMovement);
3045            if (couldChooseFirst)
3046              printf("candidate %d up %g down %g sort %g\n",iDo,choice.upMovement,choice.downMovement,sort[iDo]);
3047#if ZERO_ONE==2
3048            // branch on 0-1 first (temp)
3049            if (fabs(choice.possibleBranch->value())<1.0) {
3050              choice.upMovement *= ZERO_FAKE;
3051              choice.downMovement *= ZERO_FAKE;
3052            }
3053#endif
3054            // feasible - see which best
3055            if (!canSkip) {
3056              if (iColumn==-46) {
3057                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
3058                     upEstimate[iObject]);
3059                printf("downMove %g upMove %g value %g current pseudo %g %g\n",
3060                       choice.downMovement,choice.upMovement,choice.possibleBranch->value(),
3061                       dynamicObject->downDynamicPseudoCost(),dynamicObject->upDynamicPseudoCost());
3062              }
3063              if (model->messageHandler()->logLevel()>3) 
3064                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
3065                     upEstimate[iObject]);
3066              model->messageHandler()->message(CBC_STRONG,model->messages())
3067                << iObject << iColumn
3068                <<choice.downMovement<<choice.numIntInfeasDown
3069                <<choice.upMovement<<choice.numIntInfeasUp
3070                <<choice.possibleBranch->value()
3071                <<CoinMessageEol;
3072            }
3073            //if (!stateOfSearch)
3074            //choice.numIntInfeasDown=99999; // temp fudge
3075            if (wantMiniTree)
3076              decision->setBestCriterion(-1.0);
3077            double bestCriterion = -1.0;
3078            //double gap = saveUpper[iColumn]-saveLower[iColumn];
3079            // Give precedence to ones with gap of 1.0
3080            //assert(gap>0.0);
3081            double factor = 1.0; //changeFactor/CoinMin(gap,100.0);
3082            int betterWay = decision->betterBranch(choice.possibleBranch,
3083                                                   branch_,
3084                                                   choice.upMovement*factor,
3085                                                   choice.numIntInfeasUp ,
3086                                                   choice.downMovement*factor,
3087                                                   choice.numIntInfeasDown );
3088            if (wantMiniTree) {
3089              double criterion = decision->getBestCriterion();
3090              sort[numberMini]=-criterion;
3091              whichObject[numberMini++]=whichObject[iDo];
3092              assert (betterWay);
3093              if (criterion>bestCriterion) 
3094                bestCriterion=criterion;
3095              else
3096                betterWay=0;
3097            }
3098            if (iDo>=changeStrategy) {
3099              // make less likely
3100              changeStrategy+=numberStrong;
3101              changeFactor *= 0.9;
3102            }
3103            if (betterWay) {
3104              delete branch_;
3105              // C) create branching object
3106              branch_ = choice.possibleBranch;
3107              choice.possibleBranch=NULL;
3108              branch_->way(betterWay);
3109              if (couldChooseFirst)
3110                printf("choosing %d way %d\n",iDo,betterWay);
3111              bestChoice = choice.objectNumber;
3112              whichChoice = iDo;
3113              if (numberStrong<=1) {
3114                delete ws;
3115                ws=NULL;
3116                break;
3117              }
3118            } else {
3119              delete choice.possibleBranch;
3120              choice.possibleBranch=NULL;
3121              if (iDo>=2*numberStrong) {
3122                delete ws;
3123                ws=NULL;
3124                break;
3125              }
3126              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
3127                if (iDo-whichChoice>=numberStrong) {
3128                  delete choice.possibleBranch;
3129                  choice.possibleBranch=NULL;
3130                  break; // give up
3131                }
3132              } else {
3133                if (iDo-whichChoice>=2*numberStrong) {
3134                  delete ws;
3135                  ws=NULL;
3136                  delete choice.possibleBranch;
3137                  choice.possibleBranch=NULL;
3138                  break; // give up
3139                }
3140              }
3141            }
3142          } else {
3143            // up feasible, down infeasible
3144            anyAction=-1;
3145            worstFeasible = CoinMax(worstFeasible,choice.upMovement);
3146            //printf("Down infeasible for choice %d sequence %d\n",i,
3147            // model->object(choice.objectNumber)->columnNumber());
3148            if (!solveAll) {
3149              choice.possibleBranch->way(1);
3150              choice.possibleBranch->branch();
3151              delete choice.possibleBranch;
3152              choice.possibleBranch=NULL;
3153              delete ws;
3154              ws=NULL;
3155              break;
3156            } else {
3157              choice.fix=1;
3158              fixObject[numberToFix++]=choice;
3159              choice.possibleBranch=NULL;
3160#define FIXNOW
3161#ifdef FIXNOW
3162              double value = ceil(saveSolution[iColumn]);
3163              saveLower[iColumn]=value;
3164              solver->setColLower(iColumn,value);
3165              assert(doneHotStart);
3166              solver->unmarkHotStart();
3167              solver->markHotStart();
3168              // may be infeasible (if other way stopped on iterations)
3169              if (!solver->isProvenOptimal()) {
3170                // neither side feasible
3171                anyAction=-2;
3172                delete choice.possibleBranch;
3173                choice.possibleBranch=NULL;
3174                //printf("Both infeasible for choice %d sequence %d\n",i,
3175                // model->object(choice.objectNumber)->columnNumber());
3176                delete ws;
3177                ws=NULL;
3178                break;
3179              }
3180#endif
3181            }
3182          }
3183        } else {
3184          if(choice.downMovement<1.0e100) {
3185            // down feasible, up infeasible
3186            anyAction=-1;
3187            worstFeasible = CoinMax(worstFeasible,choice.downMovement);
3188            //printf("Up infeasible for choice %d sequence %d\n",i,
3189            // model->object(choice.objectNumber)->columnNumber());
3190            if (!solveAll) {
3191              choice.possibleBranch->way(-1);
3192              choice.possibleBranch->branch();
3193              delete choice.possibleBranch;
3194              choice.possibleBranch=NULL;
3195              delete ws;
3196              ws=NULL;
3197              break;
3198            } else {
3199              choice.fix=-1;
3200              fixObject[numberToFix++]=choice;
3201              choice.possibleBranch=NULL;
3202#ifdef FIXNOW
3203              double value = floor(saveSolution[iColumn]);
3204              saveUpper[iColumn]=value;
3205              solver->setColUpper(iColumn,value);
3206              assert(doneHotStart);
3207              solver->unmarkHotStart();
3208              solver->markHotStart();
3209              // may be infeasible (if other way stopped on iterations)
3210              if (!solver->isProvenOptimal()) {
3211                // neither side feasible
3212                anyAction=-2;
3213                delete choice.possibleBranch;
3214                choice.possibleBranch=NULL;
3215                //printf("Both infeasible for choice %d sequence %d\n",i,
3216                // model->object(choice.objectNumber)->columnNumber());
3217                delete ws;
3218                ws=NULL;
3219                break;
3220              }
3221#endif
3222            }
3223          } else {
3224            // neither side feasible
3225            anyAction=-2;
3226            delete choice.possibleBranch;
3227            choice.possibleBranch=NULL;
3228            //printf("Both infeasible for choice %d sequence %d\n",i,
3229            // model->object(choice.objectNumber)->columnNumber());
3230            delete ws;
3231            ws=NULL;
3232            break;
3233          }
3234        }
3235        // Check max time
3236        hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
3237                       model->getDblParam(CbcModel::CbcMaximumSeconds));
3238        if (hitMaxTime) {
3239          // make sure rest are fast
3240          doQuickly=true;
3241          for ( int jDo=iDo+1;jDo<numberToDo;jDo++) {
3242            int iObject = whichObject[iDo];
3243            CbcObject * object = model->modifiableObject(iObject);
3244            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3245              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3246            dynamicObject->setNumberBeforeTrust(0);
3247          }
3248          numberTest=0;
3249          distanceToCutoff=-COIN_DBL_MAX;
3250        }
3251        delete choice.possibleBranch;
3252      }
3253      double averageChange = model->sumChangeObjective()/((double) model->getNodeCount());
3254      if (depth_<10||worstFeasible>0.2*averageChange) 
3255        solveAll=false;
3256      if (model->messageHandler()->logLevel()>3||false) { 
3257        if (anyAction==-2)
3258          printf("infeasible\n");
3259        else if(anyAction==-1)
3260          if (!solveAll)
3261            printf("%d fixed\n",numberToFix);
3262          else
3263            printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n",numberToFix,bestChoice,
3264                   iDo,whichChoice,numberToDo);
3265        else
3266          printf("choosing %d iDo %d iChosenWhen %d numberToDo %d\n",bestChoice,
3267                 iDo,whichChoice,numberToDo);
3268      }
3269      if (doneHotStart) {
3270        // Delete the snapshot
3271        solver->unmarkHotStart();
3272        // back to normal
3273        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3274        // restore basis
3275        solver->setWarmStart(ws);
3276      }
3277      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
3278      // Unless infeasible we will carry on
3279      // But we could fix anyway
3280      if (numberToFix&&!hitMaxTime) {
3281        if (anyAction==-2) {
3282          // take off
3283          for (i = 0 ; i < numberToFix ; i++) {
3284            delete fixObject[i].possibleBranch;
3285          }
3286        } else {
3287          // apply and take off
3288          for (i = 0 ; i < numberToFix ; i++) {
3289#ifndef FIXNOW
3290            fixObject[i].possibleBranch->way(fixObject[i].fix) ;
3291            fixObject[i].possibleBranch->branch() ;
3292#endif
3293            delete fixObject[i].possibleBranch;
3294          }
3295          bool feasible=true;
3296#if ACTION <2
3297          if (solveAll) {
3298            // can do quick optimality check
3299            int easy=2;
3300            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3301            solver->resolve() ;
3302            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3303            feasible = solver->isProvenOptimal();
3304            if (feasible) {
3305              anyAction=0;
3306              numberMini=0;
3307              memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
3308              model->reserveCurrentSolution(saveSolution);
3309              memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
3310              memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
3311              model->setPointers(solver);
3312              // See if candidate still possible
3313              if (branch_) {
3314                const CbcObject * object = model->object(bestChoice);
3315                int preferredWay;
3316                double infeasibility = object->infeasibility(preferredWay);
3317                if (!infeasibility) {
3318                  // take out
3319                  delete branch_;
3320                  branch_=NULL;
3321                } else {
3322                  branch_->way(preferredWay);
3323                }
3324              }
3325            } else {
3326              anyAction=-2;
3327              finished=true;
3328            }
3329          }
3330#endif
3331          // If  fixed then round again
3332          if (!branch_&&anyAction!=-2) {
3333            finished=false;
3334          }
3335          // If these in then different action
3336#if ACTION == 1
3337          if (!anyAction)
3338            anyAction=-1;
3339          finished=true;
3340#endif
3341        }
3342      }
3343      delete ws;
3344    }
3345  }
3346  if (model->messageHandler()->logLevel()>2) 
3347    printf("%d strong, %d iters, %d pen, %d mark, %d fixed, action %d nnott %d nt %d, %d dq %s ns %d\n",
3348         numberStrongDone,numberStrongIterations,xPen,xMark,
3349           numberToFix,anyAction,numberNotTrusted,px[0],px[1],px[2]>0 ? "y" : "n",px[3]);
3350  // update number of strong iterations etc
3351  model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
3352                             anyAction==-2 ? 0:numberToFix,anyAction==-2);
3353  if (!newWay) {
3354  if (((model->searchStrategy()+1)%1000)==0) {
3355    if (solver->messageHandler()->logLevel()>1)
3356      printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3357             numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3358             numberNotTrusted);
3359    // decide what to do
3360    int strategy=1;
3361    if (numberUnfinished*4>numberStrongDone&&numberStrongInfeasible*10<numberStrongDone) {
3362      strategy=2;
3363      if (model->logLevel()>1)
3364        printf("going to strategy 2\n");
3365    }
3366    if (numberNodes)
3367      strategy=1;  // should only happen after hot start
3368    model->setSearchStrategy(strategy);
3369  }
3370  }
3371  //if (numberToFix&&depth_<5)
3372  //printf("%d fixed by strong at depth %d\n",numberToFix,depth_);
3373  // Set guessed solution value
3374  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
3375 
3376  // Get collection of branches if mini tree wanted
3377  if (anyAction==0&&numberMini&&numberMini>1) {
3378    // Sort
3379    CoinSort_2(sort,sort+numberMini,whichObject);
3380    delete branch_;
3381    branch_=NULL;
3382    numberMini = CoinMin(numberMini,model->sizeMiniTree());
3383    anyAction=numberMini;
3384    branches = new OsiSolverBranch[numberMini];
3385    for (int iDo=0;iDo<numberMini;iDo++) {
3386      int iObject = whichObject[iDo];
3387      CbcObject * object = model->modifiableObject(iObject);
3388      OsiSolverBranch * oneBranch = object->solverBranch();
3389      branches[iDo]=*oneBranch;
3390      delete oneBranch;
3391    }
3392  }
3393/*
3394  Cleanup, then we're finished
3395*/
3396  if (!model->branchingMethod())
3397    delete decision;
3398   
3399  delete [] fixObject;
3400  delete [] sort;
3401  delete [] whichObject;
3402  delete [] objectMark;
3403  delete [] saveLower;
3404  delete [] saveUpper;
3405  delete [] upEstimate;
3406  delete [] downEstimate;
3407# ifdef COIN_HAS_CLP
3408  if (osiclp) 
3409    osiclp->setSpecialOptions(saveClpOptions);
3410# endif
3411  // restore solution
3412  solver->setColSolution(saveSolution);
3413  model->reserveCurrentSolution(saveSolution);
3414  delete [] saveSolution;
3415  model->setStateOfSearch(saveStateOfSearch);
3416  model->setLogLevel(saveLogLevel);
3417  return anyAction;
3418}
3419int CbcNode::analyze (CbcModel *model, double * results)
3420{
3421  int i;
3422  int numberIterationsAllowed = model->numberAnalyzeIterations();
3423  OsiSolverInterface * solver = model->solver();
3424  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
3425  double cutoff =model->getCutoff();
3426  const double * lower = solver->getColLower();
3427  const double * upper = solver->getColUpper();
3428  const double * dj = solver->getReducedCost();
3429  int numberObjects = model->numberObjects();
3430  int numberColumns = model->getNumCols();
3431  // Initialize arrays
3432  int numberIntegers = model->numberIntegers();
3433  int * back = new int[numberColumns];
3434  const int * integerVariable = model->integerVariable();
3435  for (i=0;i<numberColumns;i++) 
3436    back[i]=-1;
3437  // What results is
3438  double * newLower = results;
3439  double * objLower = newLower+numberIntegers;
3440  double * newUpper = objLower+numberIntegers;
3441  double * objUpper = newUpper+numberIntegers;
3442  for (i=0;i<numberIntegers;i++) {
3443    int iColumn = integerVariable[i];
3444    back[iColumn]=i;
3445    newLower[i]=0.0;
3446    objLower[i]=-COIN_DBL_MAX;
3447    newUpper[i]=0.0;
3448    objUpper[i]=-COIN_DBL_MAX;
3449  }
3450  double * saveUpper = new double[numberColumns];
3451  double * saveLower = new double[numberColumns];
3452  int anyAction=0;
3453  // Save solution in case heuristics need good solution later
3454 
3455  double * saveSolution = new double[numberColumns];
3456  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
3457  model->reserveCurrentSolution(saveSolution);
3458  for (i=0;i<numberColumns;i++) {
3459    saveLower[i] = lower[i];
3460    saveUpper[i] = upper[i];
3461  }
3462  // Get arrays to sort
3463  double * sort = new double[numberObjects];
3464  int * whichObject = new int[numberObjects];
3465  int numberToFix=0;
3466  int numberToDo=0;
3467     
3468  // compute current state
3469  int numberObjectInfeasibilities; // just odd ones
3470  int numberIntegerInfeasibilities;
3471  model->feasibleSolution(
3472                          numberIntegerInfeasibilities,
3473                          numberObjectInfeasibilities);
3474# ifdef COIN_HAS_CLP
3475  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
3476  int saveClpOptions=0;
3477  bool fastIterations = (model->specialOptions()&8)!=0;
3478  if (osiclp&&fastIterations) {
3479    // for faster hot start
3480    saveClpOptions = osiclp->specialOptions();
3481    osiclp->setSpecialOptions(saveClpOptions|8192);
3482  }
3483# else
3484  bool fastIterations = false ;
3485# endif
3486  /*
3487    Scan for branching objects that indicate infeasibility. Choose candidates
3488    using priority as the first criteria, then integer infeasibility.
3489   
3490    The algorithm is to fill the array with a set of good candidates (by
3491    infeasibility) with priority bestPriority.  Finding a candidate with
3492    priority better (less) than bestPriority flushes the choice array. (This
3493    serves as initialization when the first candidate is found.)
3494   
3495  */
3496  numberToDo=0;
3497  for (i=0;i<numberObjects;i++) {
3498    CbcObject * object = model->modifiableObject(i);
3499    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3500      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3501    if(!dynamicObject)
3502      continue;
3503    int preferredWay;
3504    double infeasibility = object->infeasibility(preferredWay);
3505    int iColumn = dynamicObject->columnNumber();
3506    if (saveUpper[iColumn]==saveLower[iColumn])
3507      continue;
3508    if (infeasibility)
3509      sort[numberToDo]=-1.0e10-infeasibility;
3510    else
3511      sort[numberToDo]=-fabs(dj[iColumn]);
3512    whichObject[numberToDo++]=i;
3513  }
3514  // Save basis
3515  CoinWarmStart * ws = solver->getWarmStart();
3516  int saveLimit;
3517  solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
3518  int targetIterations = CoinMax(500,numberIterationsAllowed/numberObjects);
3519  if (saveLimit<targetIterations)
3520    solver->setIntParam(OsiMaxNumIterationHotStart,targetIterations); 
3521  // Mark hot start
3522  solver->markHotStart();
3523  // Sort
3524  CoinSort_2(sort,sort+numberToDo,whichObject);
3525  //double distanceToCutoff=model->getCutoff()-objectiveValue_;
3526  double * currentSolution = model->currentSolution();
3527  double integerTolerance = 
3528    model->getDblParam(CbcModel::CbcIntegerTolerance);
3529  double objMin = 1.0e50;
3530  double objMax = -1.0e50;
3531  bool needResolve=false;
3532  int iDo;
3533  for (iDo=0;iDo<numberToDo;iDo++) {
3534    CbcStrongInfo choice;
3535    int iObject = whichObject[iDo];
3536    CbcObject * object = model->modifiableObject(iObject);
3537    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3538      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3539    int iColumn = dynamicObject->columnNumber();
3540    int preferredWay;
3541    object->infeasibility(preferredWay);
3542    double value = currentSolution[iColumn];
3543    double nearest = floor(value+0.5);
3544    double lowerValue = floor(value);
3545    bool satisfied=false;
3546    if (fabs(value-nearest)<=integerTolerance||value<saveLower[iColumn]||value>saveUpper[iColumn]) {
3547      satisfied=true;
3548      double newValue;
3549      if (nearest<saveUpper[iColumn]) {
3550        newValue = nearest + 1.0001*integerTolerance;
3551        lowerValue = nearest;
3552      } else {
3553        newValue = nearest - 1.0001*integerTolerance;
3554        lowerValue = nearest-1;
3555      }
3556      currentSolution[iColumn]=newValue;
3557    }
3558    double upperValue = lowerValue+1.0;
3559    choice.possibleBranch=object->createBranch(preferredWay);
3560    currentSolution[iColumn]=value;
3561    // Save which object it was
3562    choice.objectNumber=iObject;
3563    choice.numIntInfeasUp = numberUnsatisfied_;
3564    choice.numIntInfeasDown = numberUnsatisfied_;
3565    choice.downMovement = 0.0;
3566    choice.upMovement = 0.0;
3567    choice.numItersDown = 0;
3568    choice.numItersUp = 0;
3569    choice.fix=0; // say not fixed
3570    double objectiveChange ;
3571    double newObjectiveValue=1.0e100;
3572    int j;
3573    // status is 0 finished, 1 infeasible and other
3574    int iStatus;
3575    /*
3576      Try the down direction first. (Specify the initial branching alternative as
3577      down with a call to way(-1). Each subsequent call to branch() performs the
3578      specified branch and advances the branch object state to the next branch
3579      alternative.)
3580    */
3581    choice.possibleBranch->way(-1) ;
3582    choice.possibleBranch->branch() ;
3583    if (fabs(value-lowerValue)>integerTolerance) {
3584      solver->solveFromHotStart() ;
3585      /*
3586        We now have an estimate of objective degradation that we can use for strong
3587        branching. If we're over the cutoff, the variable is monotone up.
3588        If we actually made it to optimality, check for a solution, and if we have
3589        a good one, call setBestSolution to process it. Note that this may reduce the
3590        cutoff, so we check again to see if we can declare this variable monotone.
3591      */
3592      if (solver->isProvenOptimal())
3593        iStatus=0; // optimal
3594      else if (solver->isIterationLimitReached()
3595               &&!solver->isDualObjectiveLimitReached())
3596        iStatus=2; // unknown
3597      else
3598        iStatus=1; // infeasible
3599      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3600      choice.numItersDown = solver->getIterationCount();
3601      numberIterationsAllowed -= choice.numItersDown;
3602      objectiveChange = newObjectiveValue  - objectiveValue_;
3603      if (!iStatus) {
3604        choice.finishedDown = true ;
3605        if (newObjectiveValue>=cutoff) {
3606          objectiveChange = 1.0e100; // say infeasible
3607        } else {
3608          // See if integer solution
3609          if (model->feasibleSolution(choice.numIntInfeasDown,
3610                                      choice.numObjInfeasDown)
3611              &&model->problemFeasibility()->feasible(model,-1)>=0) {
3612            model->setBestSolution(CBC_STRONGSOL,
3613                                   newObjectiveValue,
3614                                   solver->getColSolution()) ;
3615            model->setLastHeuristic(NULL);
3616            model->incrementUsed(solver->getColSolution());
3617            cutoff =model->getCutoff();
3618            if (newObjectiveValue >= cutoff)    //  *new* cutoff
3619              objectiveChange = 1.0e100 ;
3620          }
3621        }
3622      } else if (iStatus==1) {
3623        objectiveChange = 1.0e100 ;
3624      } else {
3625        // Can't say much as we did not finish
3626        choice.finishedDown = false ;
3627      }
3628      choice.downMovement = objectiveChange ;
3629    }
3630    // restore bounds
3631    for ( j=0;j<numberColumns;j++) {
3632      if (saveLower[j] != lower[j])
3633        solver->setColLower(j,saveLower[j]);
3634      if (saveUpper[j] != upper[j])
3635        solver->setColUpper(j,saveUpper[j]);
3636    }
3637    // repeat the whole exercise, forcing the variable up
3638    choice.possibleBranch->branch();
3639    if (fabs(value-upperValue)>integerTolerance) {
3640      solver->solveFromHotStart() ;
3641      /*
3642        We now have an estimate of objective degradation that we can use for strong
3643        branching. If we're over the cutoff, the variable is monotone up.
3644        If we actually made it to optimality, check for a solution, and if we have
3645        a good one, call setBestSolution to process it. Note that this may reduce the
3646        cutoff, so we check again to see if we can declare this variable monotone.
3647      */
3648      if (solver->isProvenOptimal())
3649        iStatus=0; // optimal
3650      else if (solver->isIterationLimitReached()
3651               &&!solver->isDualObjectiveLimitReached())
3652        iStatus=2; // unknown
3653      else
3654        iStatus=1; // infeasible
3655      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3656      choice.numItersUp = solver->getIterationCount();
3657      numberIterationsAllowed -= choice.numItersUp;
3658      objectiveChange = newObjectiveValue  - objectiveValue_;
3659      if (!iStatus) {
3660        choice.finishedUp = true ;
3661        if (newObjectiveValue>=cutoff) {
3662          objectiveChange = 1.0e100; // say infeasible
3663        } else {
3664          // See if integer solution
3665          if (model->feasibleSolution(choice.numIntInfeasUp,
3666                                      choice.numObjInfeasUp)
3667              &&model->problemFeasibility()->feasible(model,-1)>=0) {
3668            model->setBestSolution(CBC_STRONGSOL,
3669                                   newObjectiveValue,
3670                                   solver->getColSolution()) ;
3671            model->setLastHeuristic(NULL);
3672            model->incrementUsed(solver->getColSolution());
3673            cutoff =model->getCutoff();
3674            if (newObjectiveValue >= cutoff)    //  *new* cutoff
3675              objectiveChange = 1.0e100 ;
3676          }
3677        }
3678      } else if (iStatus==1) {
3679        objectiveChange = 1.0e100 ;
3680      } else {
3681        // Can't say much as we did not finish
3682        choice.finishedUp = false ;
3683      }
3684      choice.upMovement = objectiveChange ;
3685     
3686      // restore bounds
3687      for ( j=0;j<numberColumns;j++) {
3688        if (saveLower[j] != lower[j])
3689          solver->setColLower(j,saveLower[j]);
3690        if (saveUpper[j] != upper[j])
3691          solver->setColUpper(j,saveUpper[j]);
3692      }
3693    }
3694    // If objective goes above certain amount we can set bound
3695    int jInt = back[iColumn];
3696    newLower[jInt]=upperValue;
3697    if (choice.finishedDown||!fastIterations)
3698      objLower[jInt]=choice.downMovement+objectiveValue_;
3699    else
3700      objLower[jInt]=objectiveValue_;
3701    newUpper[jInt]=lowerValue;
3702    if (choice.finishedUp||!fastIterations)
3703      objUpper[jInt]=choice.upMovement+objectiveValue_;
3704    else
3705      objUpper[jInt]=objectiveValue_;
3706    objMin = CoinMin(CoinMin(objLower[jInt],objUpper[jInt]),objMin);
3707    /*
3708      End of evaluation for this candidate variable. Possibilities are:
3709      * Both sides below cutoff; this variable is a candidate for branching.
3710      * Both sides infeasible or above the objective cutoff: no further action
3711      here. Break from the evaluation loop and assume the node will be purged
3712      by the caller.
3713      * One side below cutoff: Install the branch (i.e., fix the variable). Break
3714      from the evaluation loop and assume the node will be reoptimised by the
3715      caller.
3716    */
3717    if (choice.upMovement<1.0e100) {
3718      if(choice.downMovement<1.0e100) {
3719        objMax = CoinMax(CoinMax(objLower[jInt],objUpper[jInt]),objMax);
3720        // In case solution coming in was odd
3721        choice.upMovement = CoinMax(0.0,choice.upMovement);
3722        choice.downMovement = CoinMax(0.0,choice.downMovement);
3723        // feasible -
3724        model->messageHandler()->message(CBC_STRONG,model->messages())
3725          << iObject << iColumn
3726          <<choice.downMovement<<choice.numIntInfeasDown
3727          <<choice.upMovement<<choice.numIntInfeasUp
3728          <<value
3729          <<CoinMessageEol;
3730      } else {
3731        // up feasible, down infeasible
3732        anyAction=-1;
3733        if (!satisfied)
3734          needResolve=true;
3735        choice.fix=1;
3736        numberToFix++;
3737        saveLower[iColumn]=upperValue;
3738        solver->setColLower(iColumn,upperValue);
3739      }
3740    } else {
3741      if(choice.downMovement<1.0e100) {
3742        // down feasible, up infeasible
3743        anyAction=-1;
3744        if (!satisfied)
3745          needResolve=true;
3746        choice.fix=-1;
3747        numberToFix++;
3748        saveUpper[iColumn]=lowerValue;
3749        solver->setColUpper(iColumn,lowerValue);
3750      } else {
3751        // neither side feasible
3752        anyAction=-2;
3753        printf("Both infeasible for choice %d sequence %d\n",i,
3754               model->object(choice.objectNumber)->columnNumber());
3755        delete ws;
3756        ws=NULL;
3757        solver->writeMps("bad");
3758        numberToFix=-1;
3759        delete choice.possibleBranch;
3760        choice.possibleBranch=NULL;
3761        break;
3762      }
3763    }
3764    delete choice.possibleBranch;
3765    if (numberIterationsAllowed<=0)
3766      break;
3767    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
3768    //     choice.downMovement,choice.upMovement,value);
3769  }
3770  printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
3771         objMin,objMax,iDo,model->numberAnalyzeIterations()-numberIterationsAllowed);
3772  model->setNumberAnalyzeIterations(numberIterationsAllowed);
3773  // Delete the snapshot
3774  solver->unmarkHotStart();
3775  // back to normal
3776  solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3777  solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
3778  // restore basis
3779  solver->setWarmStart(ws);
3780  delete ws;
3781   
3782  delete [] sort;
3783  delete [] whichObject;
3784  delete [] saveLower;
3785  delete [] saveUpper;
3786  delete [] back;
3787  // restore solution
3788  solver->setColSolution(saveSolution);
3789# ifdef COIN_HAS_CLP
3790  if (osiclp) 
3791    osiclp->setSpecialOptions(saveClpOptions);
3792# endif
3793  model->reserveCurrentSolution(saveSolution);
3794  delete [] saveSolution;
3795  if (needResolve)
3796    solver->resolve();
3797  return numberToFix;
3798}
3799
3800
3801CbcNode::CbcNode(const CbcNode & rhs) 
3802{ 
3803#ifdef CHECK_NODE
3804  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
3805#endif
3806  if (rhs.nodeInfo_)
3807    nodeInfo_ = rhs.nodeInfo_->clone();
3808  else
3809    nodeInfo_=NULL;
3810  objectiveValue_=rhs.objectiveValue_;
3811  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3812  if (rhs.branch_)
3813    branch_=rhs.branch_->clone();
3814  else
3815    branch_=NULL;
3816  depth_ = rhs.depth_;
3817  numberUnsatisfied_ = rhs.numberUnsatisfied_;
3818}
3819
3820CbcNode &
3821CbcNode::operator=(const CbcNode & rhs)
3822{
3823  if (this != &rhs) {
3824    delete nodeInfo_;
3825    if (rhs.nodeInfo_)
3826      nodeInfo_ = rhs.nodeInfo_->clone();
3827    else
3828      nodeInfo_ = NULL;
3829    objectiveValue_=rhs.objectiveValue_;
3830    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3831    if (rhs.branch_)
3832      branch_=rhs.branch_->clone();
3833    else
3834      branch_=NULL,
3835    depth_ = rhs.depth_;
3836    numberUnsatisfied_ = rhs.numberUnsatisfied_;
3837  }
3838  return *this;
3839}
3840
3841
3842CbcNode::~CbcNode ()
3843{
3844#ifdef CHECK_NODE
3845  if (nodeInfo_) 
3846    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
3847         this,nodeInfo_,nodeInfo_->numberPointingToThis());
3848  else
3849    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
3850         this,nodeInfo_);
3851#endif
3852  if (nodeInfo_) {
3853    nodeInfo_->nullOwner();
3854    int numberToDelete=nodeInfo_->numberBranchesLeft();
3855    //    CbcNodeInfo * parent = nodeInfo_->parent();
3856    if (nodeInfo_->decrement(numberToDelete)==0) {
3857      delete nodeInfo_;
3858    } else {
3859      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent);
3860      // anyway decrement parent
3861      //if (parent)
3862      ///parent->decrement(1);
3863    }
3864  }
3865  delete branch_;
3866}
3867// Decrement  active cut counts
3868void 
3869CbcNode::decrementCuts(int change)
3870{
3871  if(nodeInfo_) {
3872    nodeInfo_->decrementCuts(change);
3873  }
3874}
3875void 
3876CbcNode::decrementParentCuts(int change)
3877{
3878  if(nodeInfo_) {
3879    nodeInfo_->decrementParentCuts(change);
3880  }
3881}
3882
3883/*
3884  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
3885  in the attached nodeInfo_.
3886*/
3887void
3888CbcNode::initializeInfo()
3889{
3890  assert(nodeInfo_ && branch_) ;
3891  nodeInfo_->initializeInfo(branch_->numberBranches());
3892}
3893// Nulls out node info
3894void 
3895CbcNode::nullNodeInfo()
3896{
3897  nodeInfo_=NULL;
3898}
3899
3900int
3901CbcNode::branch()
3902{
3903  double changeInGuessed=branch_->branch(true);
3904  guessedObjectiveValue_+= changeInGuessed;
3905  //#define PRINTIT
3906#ifdef PRINTIT
3907  int numberLeft = nodeInfo_->numberBranchesLeft();
3908  CbcNodeInfo * parent = nodeInfo_->parent();
3909  int parentNodeNumber = -1;
3910  //CbcBranchingObject * object1 = branch_->object_;
3911  //CbcObject * object = object1->
3912  //int sequence = object->columnNumber);
3913  int id=-1;
3914  double value=0.0;
3915  if (branch_) {
3916    id = branch_->variable();
3917    value = branch_->value();
3918  }
3919  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
3920  if (parent)
3921    parentNodeNumber = parent->nodeNumber();
3922  printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
3923         nodeInfo_->nodeNumber(),(numberLeft==2) ? "leftBranch" : "rightBranch",
3924         way(),depth_,parentNodeNumber);
3925#endif
3926  return nodeInfo_->branchedOn();
3927}
Note: See TracBrowser for help on using the repository browser.