source: trunk/Cbc/src/CbcNode.cpp @ 789

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

trying to move ampl stuff away

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