source: trunk/CbcNode.cpp @ 202

Last change on this file since 202 was 202, checked in by forrest, 16 years ago

try and improve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 98.4 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#include <string>
8//#define CBC_DEBUG 1
9//#define CHECK_CUT_COUNTS
10//#define CHECK_NODE
11#include <cassert>
12#include <cfloat>
13#define CUTS
14#include "OsiSolverInterface.hpp"
15#include "CoinWarmStartBasis.hpp"
16#include "CoinTime.hpp"
17#include "CbcModel.hpp"
18#include "CbcNode.hpp"
19#include "CbcBranchActual.hpp"
20#include "CbcBranchDynamic.hpp"
21#include "OsiRowCut.hpp"
22#include "OsiRowCutDebugger.hpp"
23#include "OsiCuts.hpp"
24#include "CbcCountRowCut.hpp"
25#include "CbcFeasibilityBase.hpp"
26#include "CbcMessage.hpp"
27#include "OsiClpSolverInterface.hpp"
28#include "ClpSimplexOther.hpp"
29using namespace std;
30#include "CglCutGenerator.hpp"
31// Default Constructor
32CbcNodeInfo::CbcNodeInfo ()
33  :
34  numberPointingToThis_(0),
35  parent_(NULL),
36  owner_(NULL),
37  numberCuts_(0),
38  nodeNumber_(0),
39  cuts_(NULL),
40  numberRows_(0),
41  numberBranchesLeft_(0)
42{
43#ifdef CHECK_NODE
44  printf("CbcNodeInfo %x Constructor\n",this);
45#endif
46}
47// Constructor given parent
48CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent)
49  :
50  numberPointingToThis_(2),
51  parent_(parent),
52  owner_(NULL),
53  numberCuts_(0),
54  nodeNumber_(0),
55  cuts_(NULL),
56  numberRows_(0),
57  numberBranchesLeft_(2)
58{
59#ifdef CHECK_NODE
60  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
61#endif
62  if (parent_) {
63    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
64  }
65}
66// Copy Constructor
67CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs)
68  :
69  numberPointingToThis_(rhs.numberPointingToThis_),
70  parent_(rhs.parent_),
71  owner_(rhs.owner_),
72  numberCuts_(rhs.numberCuts_),
73  nodeNumber_(rhs.nodeNumber_),
74  cuts_(NULL),
75  numberRows_(rhs.numberRows_),
76  numberBranchesLeft_(rhs.numberBranchesLeft_)
77{
78#ifdef CHECK_NODE
79  printf("CbcNodeInfo %x Copy constructor\n",this);
80#endif
81  if (numberCuts_) {
82    cuts_ = new CbcCountRowCut * [numberCuts_];
83    int n=0;
84    for (int i=0;i<numberCuts_;i++) {
85      CbcCountRowCut * thisCut = rhs.cuts_[i];
86      if (thisCut) {
87        // I think this is correct - new one should take priority
88        thisCut->setInfo(this,n);
89        thisCut->increment(numberBranchesLeft_); 
90        cuts_[n++] = thisCut;
91      }
92    }
93    numberCuts_=n;
94  }
95}
96// Constructor given parent and owner
97CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent, CbcNode * owner)
98  :
99  numberPointingToThis_(2),
100  parent_(parent),
101  owner_(owner),
102  numberCuts_(0),
103  nodeNumber_(0),
104  cuts_(NULL),
105  numberRows_(0),
106  numberBranchesLeft_(2)
107{
108#ifdef CHECK_NODE
109  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
110#endif
111  if (parent_) {
112    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
113  }
114}
115
116/**
117  Take care to detach from the owning CbcNode and decrement the reference
118  count in the parent.  If this is the last nodeInfo object pointing to the
119  parent, make a recursive call to delete the parent.
120*/
121CbcNodeInfo::~CbcNodeInfo()
122{
123#ifdef CHECK_NODE
124  printf("CbcNodeInfo %x Destructor parent %x\n",this,parent_);
125#endif
126
127  assert(!numberPointingToThis_);
128  // But they may be some left (max nodes?)
129  for (int i=0;i<numberCuts_;i++) 
130    delete cuts_[i];
131
132  delete [] cuts_;
133  if (owner_) 
134    owner_->nullNodeInfo();
135  if (parent_) {
136    int numberLinks = parent_->decrement();
137    if (!numberLinks) delete parent_;
138  }
139}
140
141
142//#define ALLCUTS
143void
144CbcNodeInfo::decrementCuts(int change)
145{
146  int i;
147  // get rid of all remaining if negative
148  int changeThis;
149  if (change<0)
150    changeThis = numberBranchesLeft_;
151  else
152    changeThis = change;
153 // decrement cut counts
154  for (i=0;i<numberCuts_;i++) {
155    if (cuts_[i]) {
156      int number = cuts_[i]->decrement(changeThis);
157      if (!number) {
158        //printf("info %x del cut %d %x\n",this,i,cuts_[i]);
159        delete cuts_[i];
160        cuts_[i]=NULL;
161      }
162    }
163  }
164}
165void
166CbcNodeInfo::decrementParentCuts(int change)
167{
168  if (parent_) {
169    // get rid of all remaining if negative
170    int changeThis;
171    if (change<0)
172      changeThis = numberBranchesLeft_;
173    else
174      changeThis = change;
175    int i;
176    // Get over-estimate of space needed for basis
177    CoinWarmStartBasis dummy;
178    dummy.setSize(0,numberRows_+numberCuts_);
179    buildRowBasis(dummy);
180    /* everything is zero (i.e. free) so we can use to see
181       if latest basis */
182    CbcNodeInfo * thisInfo = parent_;
183    while (thisInfo) 
184      thisInfo = thisInfo->buildRowBasis(dummy);
185    // decrement cut counts
186    thisInfo = parent_;
187    int numberRows=numberRows_;
188    while (thisInfo) {
189      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
190        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
191#ifdef ALLCUTS
192        status = CoinWarmStartBasis::isFree;
193#endif
194        if (thisInfo->cuts_[i]) {
195          int number=1;
196          if (status!=CoinWarmStartBasis::basic) {
197            // tight - drop 1 or 2
198            if (change<0)
199              number = thisInfo->cuts_[i]->decrement(changeThis);
200            else
201              number = thisInfo->cuts_[i]->decrement(change);
202          }
203          if (!number) {
204            delete thisInfo->cuts_[i];
205            thisInfo->cuts_[i]=NULL;
206          }
207        }
208      }
209      thisInfo = thisInfo->parent_;
210    }
211  }
212}
213
214void
215CbcNodeInfo::incrementParentCuts(int change)
216{
217  if (parent_) {
218    int i;
219    // Get over-estimate of space needed for basis
220    CoinWarmStartBasis dummy;
221    dummy.setSize(0,numberRows_+numberCuts_);
222    /* everything is zero (i.e. free) so we can use to see
223       if latest basis */
224    buildRowBasis(dummy);
225    CbcNodeInfo * thisInfo = parent_;
226    while (thisInfo) 
227      thisInfo = thisInfo->buildRowBasis(dummy);
228    // increment cut counts
229    thisInfo = parent_;
230    int numberRows=numberRows_;
231    while (thisInfo) {
232      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
233        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
234#ifdef ALLCUTS
235        status = CoinWarmStartBasis::isFree;
236#endif
237        if (thisInfo->cuts_[i]&&status!=CoinWarmStartBasis::basic) {
238          thisInfo->cuts_[i]->increment(change);
239        }
240      }
241      thisInfo = thisInfo->parent_;
242    }
243  }
244}
245
246/*
247  Append cuts to the cuts_ array in a nodeInfo. The initial reference count
248  is set to numberToBranchOn, which will normally be the number of arms
249  defined for the CbcBranchingObject attached to the CbcNode that owns this
250  CbcNodeInfo.
251*/
252void
253CbcNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn,
254                      int * whichGenerator)
255{
256  int numberCuts = cuts.sizeRowCuts();
257  if (numberCuts) {
258    int i;
259    if (!numberCuts_) {
260      cuts_ = new CbcCountRowCut * [numberCuts];
261    } else {
262      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
263      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
264      delete [] cuts_;
265      cuts_ = temp;
266    }
267    for (i=0;i<numberCuts;i++) {
268      CbcCountRowCut * thisCut = new CbcCountRowCut(*cuts.rowCutPtr(i),
269                                                    this,numberCuts_);
270      thisCut->increment(numberToBranchOn); 
271      cuts_[numberCuts_++] = thisCut;
272#ifdef CBC_DEBUG
273      int n=thisCut->row().getNumElements();
274#if CBC_DEBUG>1
275      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
276             thisCut->ub());
277#endif
278      int j;
279#if CBC_DEBUG>1
280      const int * index = thisCut->row().getIndices();
281#endif
282      const double * element = thisCut->row().getElements();
283      for (j=0;j<n;j++) {
284#if CBC_DEBUG>1
285        printf(" (%d,%g)",index[j],element[j]);
286#endif
287        assert(fabs(element[j])>1.00e-12);
288      }
289      printf("\n");
290#endif
291    }
292  }
293}
294
295void
296CbcNodeInfo::addCuts(int numberCuts, CbcCountRowCut ** cut, 
297                     int numberToBranchOn)
298{
299  if (numberCuts) {
300    int i;
301    if (!numberCuts_) {
302      cuts_ = new CbcCountRowCut * [numberCuts];
303    } else {
304      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
305      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
306      delete [] cuts_;
307      cuts_ = temp;
308    }
309    for (i=0;i<numberCuts;i++) {
310      CbcCountRowCut * thisCut = cut[i];
311      thisCut->setInfo(this,numberCuts_);
312      //printf("info %x cut %d %x\n",this,i,thisCut);
313      thisCut->increment(numberToBranchOn); 
314      cuts_[numberCuts_++] = thisCut;
315#ifdef CBC_DEBUG
316      int n=thisCut->row().getNumElements();
317#if CBC_DEBUG>1
318      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
319             thisCut->ub());
320#endif
321      int j;
322#if CBC_DEBUG>1
323      const int * index = thisCut->row().getIndices();
324#endif
325      const double * element = thisCut->row().getElements();
326      for (j=0;j<n;j++) {
327#if CBC_DEBUG>1
328        printf(" (%d,%g)",index[j],element[j]);
329#endif
330        assert(fabs(element[j])>1.00e-12);
331      }
332      printf("\n");
333#endif
334    }
335  }
336}
337
338// delete cuts
339void
340CbcNodeInfo::deleteCuts(int numberToDelete, CbcCountRowCut ** cuts)
341{
342  int i;
343  int j;
344  int last=-1;
345  for (i=0;i<numberToDelete;i++) {
346    CbcCountRowCut * next = cuts[i];
347    for (j=last+1;j<numberCuts_;j++) {
348      if (next==cuts_[j])
349        break;
350    }
351    if (j==numberCuts_) {
352      // start from beginning
353      for (j=0;j<last;j++) {
354        if (next==cuts_[j])
355          break;
356      }
357      assert(j<last);
358    }
359    last=j;
360    int number = cuts_[j]->decrement();
361    if (!number)
362      delete cuts_[j];
363    cuts_[j]=NULL;
364  }
365  j=0;
366  for (i=0;i<numberCuts_;i++) {
367    if (cuts_[i])
368      cuts_[j++]=cuts_[i];
369  }
370  numberCuts_ = j;
371}
372
373// delete cuts
374void
375CbcNodeInfo::deleteCuts(int numberToDelete, int * which)
376{
377  int i;
378  for (i=0;i<numberToDelete;i++) {
379    int iCut=which[i];
380    int number = cuts_[iCut]->decrement();
381    if (!number)
382      delete cuts_[iCut];
383    cuts_[iCut]=NULL;
384  }
385  int j=0;
386  for (i=0;i<numberCuts_;i++) {
387    if (cuts_[i])
388      cuts_[j++]=cuts_[i];
389  }
390  numberCuts_ = j;
391}
392
393// Really delete a cut
394void 
395CbcNodeInfo::deleteCut(int whichOne)
396{
397  assert(whichOne<numberCuts_);
398  cuts_[whichOne]=NULL;
399}
400
401CbcFullNodeInfo::CbcFullNodeInfo() :
402  CbcNodeInfo(),
403  basis_(),
404  numberIntegers_(0),
405  lower_(NULL),
406  upper_(NULL)
407{
408}
409CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model,
410                                 int numberRowsAtContinuous) :
411  CbcNodeInfo()
412{
413  OsiSolverInterface * solver = model->solver();
414  numberRows_ = numberRowsAtContinuous;
415  numberIntegers_ = model->numberIntegers();
416  int numberColumns = model->getNumCols();
417  lower_ = new double [numberColumns];
418  upper_ = new double [numberColumns];
419  const double * lower = solver->getColLower();
420  const double * upper = solver->getColUpper();
421  int i;
422
423  for (i=0;i<numberColumns;i++) {
424    lower_[i]=lower[i];
425    upper_[i]=upper[i];
426  }
427
428  basis_ =  dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
429}
430
431CbcFullNodeInfo::CbcFullNodeInfo(const CbcFullNodeInfo & rhs) :
432  CbcNodeInfo(rhs)
433{ 
434  basis_= dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
435  numberIntegers_=rhs.numberIntegers_;
436  lower_=NULL;
437  upper_=NULL;
438  if (rhs.lower_!=NULL) {
439    int numberColumns = basis_->getNumStructural();
440    lower_ = new double [numberColumns];
441    upper_ = new double [numberColumns];
442    assert (upper_!=NULL);
443    memcpy(lower_,rhs.lower_,numberColumns*sizeof(double));
444    memcpy(upper_,rhs.upper_,numberColumns*sizeof(double));
445  }
446}
447
448CbcNodeInfo * 
449CbcFullNodeInfo::clone() const
450{ 
451  return (new CbcFullNodeInfo(*this));
452}
453
454CbcFullNodeInfo::~CbcFullNodeInfo ()
455{
456  delete basis_ ;
457  delete [] lower_;
458  delete [] upper_;
459}
460
461/*
462   The basis supplied as a parameter is deleted and replaced with a new basis
463   appropriate for the node, and lower and upper bounds on variables are
464   reset according to the stored bounds arrays. Any cuts associated with this
465   node are added to the list in addCuts, but not actually added to the
466   constraint system in the model.
467
468   Why pass in a basis at all? The short answer is ``We need the parameter to
469   pass out a basis, so might as well use it to pass in the size.''
470   
471   A longer answer is that in practice we take a memory allocation hit up in
472   addCuts1 (the only place applyToModel is called) when we setSize() the
473   basis that's passed in. It's immediately tossed here in favour of a clone
474   of the basis attached to this nodeInfo. This can probably be fixed, given
475   a bit of thought.
476*/
477
478void CbcFullNodeInfo::applyToModel (CbcModel *model,
479                                    CoinWarmStartBasis *&basis,
480                                    CbcCountRowCut **addCuts,
481                                    int &currentNumberCuts) const 
482
483{ OsiSolverInterface *solver = model->solver() ;
484
485  // branch - do bounds
486  int i;
487  solver->setColLower(lower_);
488  solver->setColUpper(upper_);
489  int numberColumns = model->getNumCols();
490  // move basis - but make sure size stays
491  int numberRows = basis->getNumArtificial();
492  delete basis ;
493  basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
494  basis->resize(numberRows,numberColumns);
495  for (i=0;i<numberCuts_;i++) 
496    addCuts[currentNumberCuts+i]= cuts_[i];
497  currentNumberCuts += numberCuts_;
498  assert(!parent_);
499  return ;
500}
501
502/* Builds up row basis backwards (until original model).
503   Returns NULL or previous one to apply .
504   Depends on Free being 0 and impossible for cuts
505*/
506CbcNodeInfo * 
507CbcFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
508{
509  const unsigned int * saved = 
510    (const unsigned int *) basis_->getArtificialStatus();
511  unsigned int * now = 
512    (unsigned int *) basis.getArtificialStatus();
513  int number=basis_->getNumArtificial()>>4;;
514  int i;
515  for (i=0;i<number;i++) { 
516    if (!now[i])
517      now[i] = saved[i];
518  }
519  return NULL;
520}
521
522// Default constructor
523CbcPartialNodeInfo::CbcPartialNodeInfo()
524
525  : CbcNodeInfo(),
526    basisDiff_(NULL),
527    variables_(NULL),
528    newBounds_(NULL),
529    numberChangedBounds_(0)
530
531{ /* this space intentionally left blank */ }
532
533// Constructor from current state
534CbcPartialNodeInfo::CbcPartialNodeInfo (CbcNodeInfo *parent, CbcNode *owner,
535                                        int numberChangedBounds,
536                                        const int *variables,
537                                        const double *boundChanges,
538                                        const CoinWarmStartDiff *basisDiff)
539 : CbcNodeInfo(parent,owner)
540{
541  basisDiff_ = basisDiff->clone() ;
542
543  numberChangedBounds_ = numberChangedBounds;
544  variables_ = new int [numberChangedBounds_];
545  newBounds_ = new double [numberChangedBounds_];
546
547  int i ;
548  for (i=0;i<numberChangedBounds_;i++) {
549    variables_[i]=variables[i];
550    newBounds_[i]=boundChanges[i];
551  }
552}
553
554CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs)
555
556  : CbcNodeInfo(rhs.parent_)
557
558{ basisDiff_ = rhs.basisDiff_->clone() ;
559
560  numberChangedBounds_ = rhs.numberChangedBounds_;
561  variables_ = new int [numberChangedBounds_];
562  newBounds_ = new double [numberChangedBounds_];
563
564  int i ;
565  for (i=0;i<numberChangedBounds_;i++) {
566    variables_[i]=rhs.variables_[i];
567    newBounds_[i]=rhs.newBounds_[i];
568  }
569}
570
571CbcNodeInfo * 
572CbcPartialNodeInfo::clone() const
573{ 
574  return (new CbcPartialNodeInfo(*this));
575}
576
577
578CbcPartialNodeInfo::~CbcPartialNodeInfo ()
579{
580  delete basisDiff_ ;
581  delete [] variables_;
582  delete [] newBounds_;
583}
584
585
586/**
587   The basis supplied as a parameter is incrementally modified, and lower and
588   upper bounds on variables in the model are incrementally modified. Any
589   cuts associated with this node are added to the list in addCuts.
590*/
591
592void CbcPartialNodeInfo::applyToModel (CbcModel *model,
593                                       CoinWarmStartBasis *&basis,
594                                       CbcCountRowCut **addCuts,
595                                       int &currentNumberCuts) const 
596
597{ OsiSolverInterface *solver = model->solver();
598
599  basis->applyDiff(basisDiff_) ;
600
601  // branch - do bounds
602  int i;
603  for (i=0;i<numberChangedBounds_;i++) {
604    int variable = variables_[i];
605    if ((variable&0x80000000)==0) {
606      // lower bound changing
607      solver->setColLower(variable,newBounds_[i]);
608    } else {
609      // upper bound changing
610      solver->setColUpper(variable&0x7fffffff,newBounds_[i]);
611    }
612  }
613  for (i=0;i<numberCuts_;i++) 
614    addCuts[currentNumberCuts+i]= cuts_[i];
615  currentNumberCuts += numberCuts_;
616  return ;
617}
618
619/* Builds up row basis backwards (until original model).
620   Returns NULL or previous one to apply .
621   Depends on Free being 0 and impossible for cuts
622*/
623
624CbcNodeInfo * 
625CbcPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
626
627{ basis.applyDiff(basisDiff_) ;
628
629  return parent_ ; }
630
631
632CbcNode::CbcNode() :
633  nodeInfo_(NULL),
634  objectiveValue_(1.0e100),
635  guessedObjectiveValue_(1.0e100),
636  branch_(NULL),
637  depth_(-1),
638  numberUnsatisfied_(0)
639{
640#ifdef CHECK_NODE
641  printf("CbcNode %x Constructor\n",this);
642#endif
643}
644
645CbcNode::CbcNode(CbcModel * model,
646                 CbcNode * lastNode) :
647  nodeInfo_(NULL),
648  objectiveValue_(1.0e100),
649  guessedObjectiveValue_(1.0e100),
650  branch_(NULL),
651  depth_(-1),
652  numberUnsatisfied_(0)
653{
654#ifdef CHECK_NODE
655  printf("CbcNode %x Constructor from model\n",this);
656#endif
657  OsiSolverInterface * solver = model->solver();
658  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
659
660  if (lastNode)
661    lastNode->nodeInfo_->increment();
662}
663
664
665void
666CbcNode::createInfo (CbcModel *model,
667                     CbcNode *lastNode,
668                     const CoinWarmStartBasis *lastws,
669                     const double *lastLower, const double *lastUpper,
670                     int numberOldActiveCuts,int numberNewCuts)
671{ OsiSolverInterface * solver = model->solver();
672/*
673  The root --- no parent. Create full basis and bounds information.
674*/
675  if (!lastNode)
676  { nodeInfo_=new CbcFullNodeInfo(model,solver->getNumRows()); }
677/*
678  Not the root. Create an edit from the parent's basis & bound information.
679  This is not quite as straightforward as it seems. We need to reintroduce
680  cuts we may have dropped out of the basis, in the correct position, because
681  this whole process is strictly positional. Start by grabbing the current
682  basis.
683*/
684  else
685  { const CoinWarmStartBasis* ws =
686      dynamic_cast<const CoinWarmStartBasis*>(solver->getWarmStart());
687    assert(ws!=NULL); // make sure not volume
688    //int numberArtificials = lastws->getNumArtificial();
689    int numberColumns = solver->getNumCols();
690   
691    const double * lower = solver->getColLower();
692    const double * upper = solver->getColUpper();
693
694    int i;
695/*
696  Create a clone and resize it to hold all the structural constraints, plus
697  all the cuts: old cuts, both active and inactive (currentNumberCuts), and
698  new cuts (numberNewCuts).
699
700  TODO: You'd think that the set of constraints (logicals) in the expanded
701        basis should match the set represented in lastws. At least, that's
702        what I thought. But at the point I first looked hard at this bit of
703        code, it turned out that lastws was the stripped basis produced at
704        the end of addCuts(), rather than the raw basis handed back by
705        addCuts1(). The expanded basis here is equivalent to the raw basis of
706        addCuts1(). I said ``whoa, that's not good, I must have introduced a
707        bug'' and went back to John's code to see where I'd gone wrong.
708        And discovered the same `error' in his code.
709
710        After a bit of thought, my conclusion is that correctness is not
711        affected by whether lastws is the stripped or raw basis. The diffs
712        have no semantics --- just a set of changes that need to be made
713        to convert lastws into expanded. I think the only effect is that we
714        store a lot more diffs (everything in expanded that's not covered by
715        the stripped basis). But I need to give this more thought. There
716        may well be some subtle error cases.
717
718        In the mean time, I've twiddled addCuts() to set lastws to the raw
719        basis. Makes me (Lou) less nervous to compare apples to apples.
720*/
721    CoinWarmStartBasis *expanded = 
722      dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
723    int numberRowsAtContinuous = model->numberRowsAtContinuous();
724    int iFull = numberRowsAtContinuous+model->currentNumberCuts()+
725      numberNewCuts;
726    //int numberArtificialsNow = iFull;
727    //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
728    //printf("l %d full %d\n",maxBasisLength,iFull);
729    expanded->resize(iFull,numberColumns);
730#ifdef FULL_DEBUG
731    printf("Before expansion: orig %d, old %d, new %d, current %d\n",
732           numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts,
733           model->currentNumberCuts()) ;
734    ws->print();
735#endif
736/*
737  Now fill in the expanded basis. Any indices beyond nPartial must
738  be cuts created while processing this node --- they can be copied directly
739  into the expanded basis. From nPartial down, pull the status of active cuts
740  from ws, interleaving with a B entry for the deactivated (loose) cuts.
741*/
742    int numberDropped = model->currentNumberCuts()-numberOldActiveCuts;
743    int iCompact=iFull-numberDropped;
744    CbcCountRowCut ** cut = model->addedCuts();
745    int nPartial = model->currentNumberCuts()+numberRowsAtContinuous;
746    iFull--;
747    for (;iFull>=nPartial;iFull--) {
748      CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
749      //assert (status != CoinWarmStartBasis::basic); // may be permanent cut
750      expanded->setArtifStatus(iFull,status);
751    }
752    for (;iFull>=numberRowsAtContinuous;iFull--) {
753      if (cut[iFull-numberRowsAtContinuous]) {
754        CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
755        // If no cut generator being used then we may have basic variables
756        //if (model->getMaximumCutPasses()&&
757        //  status == CoinWarmStartBasis::basic)
758        //printf("cut basic\n");
759        expanded->setArtifStatus(iFull,status);
760      } else {
761        expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic);
762      }
763    }
764#ifdef FULL_DEBUG
765    printf("Expanded basis\n");
766    expanded->print() ;
767    printf("Diffing against\n") ;
768    lastws->print() ;
769#endif   
770/*
771  Now that we have two bases in proper positional correspondence, creating
772  the actual diff is dead easy.
773*/
774
775    CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
776/*
777  Diff the bound vectors. It's assumed the number of structural variables is
778  not changing. Assuming that branching objects all involve integer variables,
779  we should see at least one bound change as a consequence of processing this
780  subproblem. Different types of branching objects could break this assertion.
781  Not true at all - we have not applied current branch - JJF.
782*/
783    double *boundChanges = new double [2*numberColumns] ;
784    int *variables = new int [2*numberColumns] ;
785    int numberChangedBounds=0;
786    for (i=0;i<numberColumns;i++) {
787      if (lower[i]!=lastLower[i]) {
788        variables[numberChangedBounds]=i;
789        boundChanges[numberChangedBounds++]=lower[i];
790      }
791      if (upper[i]!=lastUpper[i]) {
792        variables[numberChangedBounds]=i|0x80000000;
793        boundChanges[numberChangedBounds++]=upper[i];
794      }
795#ifdef CBC_DEBUG
796      if (lower[i]!=lastLower[i])
797        printf("lower on %d changed from %g to %g\n",
798               i,lastLower[i],lower[i]);
799      if (upper[i]!=lastUpper[i])
800        printf("upper on %d changed from %g to %g\n",
801               i,lastUpper[i],upper[i]);
802#endif
803    }
804#ifdef CBC_DEBUG
805    printf("%d changed bounds\n",numberChangedBounds) ;
806#endif
807    //if (lastNode->branchingObject()->boundBranch())
808    //assert (numberChangedBounds);
809/*
810  Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
811  return.
812*/
813    nodeInfo_ =
814      new CbcPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds,
815                             variables,boundChanges,basisDiff) ;
816    delete basisDiff ;
817    delete [] boundChanges;
818    delete [] variables;
819    delete expanded ;
820    delete ws;
821  }
822  // Set node number
823  nodeInfo_->setNodeNumber(model->getNodeCount2());
824}
825
826/*
827  The routine scans through the object list of the model looking for objects
828  that indicate infeasibility. It tests each object using strong branching
829  and selects the one with the least objective degradation.  A corresponding
830  branching object is left attached to lastNode.
831
832  If strong branching is disabled, a candidate object is chosen essentially
833  at random (whatever object ends up in pos'n 0 of the candidate array).
834
835  If a branching candidate is found to be monotone, bounds are set to fix the
836  variable and the routine immediately returns (the caller is expected to
837  reoptimize).
838
839  If a branching candidate is found to result in infeasibility in both
840  directions, the routine immediately returns an indication of infeasibility.
841
842  Returns:  0   both branch directions are feasible
843           -1   branching variable is monotone
844           -2   infeasible
845
846  Original comments:
847    Here could go cuts etc etc
848    For now just fix on objective from strong branching.
849*/
850
851int CbcNode::chooseBranch (CbcModel *model, CbcNode *lastNode,int numberPassesLeft)
852
853{ if (lastNode)
854    depth_ = lastNode->depth_+1;
855  else
856    depth_ = 0;
857  delete branch_;
858  branch_=NULL;
859  OsiSolverInterface * solver = model->solver();
860  double saveObjectiveValue = solver->getObjValue();
861  double objectiveValue = solver->getObjSense()*saveObjectiveValue;
862  const double * lower = solver->getColLower();
863  const double * upper = solver->getColUpper();
864  // See what user thinks
865  int anyAction=model->problemFeasibility()->feasible(model,0);
866  if (anyAction) {
867    // will return -2 if infeasible , 0 if treat as integer
868    return anyAction-1;
869  }
870  double integerTolerance = 
871    model->getDblParam(CbcModel::CbcIntegerTolerance);
872  int i;
873  bool beforeSolution = model->getSolutionCount()==0;
874  int numberStrong=model->numberStrong();
875  int saveNumberStrong=numberStrong;
876  int numberObjects = model->numberObjects();
877  bool checkFeasibility = numberObjects>model->numberIntegers();
878  int maximumStrong = CoinMax(CoinMin(model->numberStrong(),numberObjects),1);
879  int numberColumns = model->getNumCols();
880  double * saveUpper = new double[numberColumns];
881  double * saveLower = new double[numberColumns];
882
883  // Save solution in case heuristics need good solution later
884 
885  double * saveSolution = new double[numberColumns];
886  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
887  model->reserveCurrentSolution(saveSolution);
888  /*
889    Get a branching decision object. Use the default decision criteria unless
890    the user has loaded a decision method into the model.
891  */
892  CbcBranchDecision *decision = model->branchingMethod();
893  if (!decision)
894    decision = new CbcBranchDefaultDecision();
895
896  CbcStrongInfo * choice = new CbcStrongInfo[maximumStrong];
897  for (i=0;i<numberColumns;i++) {
898    saveLower[i] = lower[i];
899    saveUpper[i] = upper[i];
900  }
901  // May go round twice if strong branching fixes all local candidates
902  bool finished=false;
903  double estimatedDegradation=0.0; 
904  while(!finished) {
905    finished=true;
906    // Some objects may compute an estimate of best solution from here
907    estimatedDegradation=0.0; 
908    //int numberIntegerInfeasibilities=0; // without odd ones
909   
910    // We may go round this loop twice (only if we think we have solution)
911    for (int iPass=0;iPass<2;iPass++) {
912     
913      // compute current state
914      //int numberObjectInfeasibilities; // just odd ones
915      //model->feasibleSolution(
916      //                      numberIntegerInfeasibilities,
917      //                      numberObjectInfeasibilities);
918      // If forcePriority > 0 then we want best solution
919      const double * bestSolution = NULL;
920      int hotstartStrategy=model->getHotstartStrategy();
921      if (hotstartStrategy>0) {
922        bestSolution = model->bestSolution();
923      }
924     
925      // Some objects may compute an estimate of best solution from here
926      estimatedDegradation=0.0; 
927      numberUnsatisfied_ = 0;
928      int bestPriority=INT_MAX;
929      /*
930        Scan for branching objects that indicate infeasibility. Choose the best
931        maximumStrong candidates, using priority as the first criteria, then
932        integer infeasibility.
933       
934        The algorithm is to fill the choice array with a set of good candidates (by
935        infeasibility) with priority bestPriority.  Finding a candidate with
936        priority better (less) than bestPriority flushes the choice array. (This
937        serves as initialization when the first candidate is found.)
938       
939        A new candidate is added to choices only if its infeasibility exceeds the
940        current max infeasibility (mostAway). When a candidate is added, it
941        replaces the candidate with the smallest infeasibility (tracked by
942        iSmallest).
943      */
944      int iSmallest = 0;
945      double mostAway = 1.0e-100;
946      for (i = 0 ; i < maximumStrong ; i++)
947        choice[i].possibleBranch = NULL ;
948      numberStrong=0;
949      for (i=0;i<numberObjects;i++) {
950        CbcObject * object = model->modifiableObject(i);
951        int preferredWay;
952        double infeasibility = object->infeasibility(preferredWay);
953        int priorityLevel = object->priority();
954        if (bestSolution) {
955          // we are doing hot start
956          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
957          if (thisOne) {
958            int iColumn = thisOne->modelSequence();
959            if (saveUpper[iColumn]>saveLower[iColumn]) {
960              double value = saveSolution[iColumn];
961              double targetValue = bestSolution[iColumn];
962              //double originalLower = thisOne->originalLower();
963              //double originalUpper = thisOne->originalUpper();
964              // switch off if not possible
965              if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
966                /* priority outranks rest always if hotstartStrategy >1
967                   otherwise can be downgraded if at correct level.
968                   Infeasibility may be increased by targetValue to choose 1.0 values first.
969                */
970                if (fabs(value-targetValue)>integerTolerance) {
971                  if (value>targetValue) {
972                    infeasibility += value;
973                    preferredWay=-1;
974                  } else {
975                    infeasibility += targetValue;
976                    preferredWay=1;
977                  }
978                } else if (hotstartStrategy>1) {
979                  if (targetValue==saveLower[iColumn]) {
980                    infeasibility += integerTolerance+1.0e-12;
981                    preferredWay=-1;
982                  } else if (targetValue==saveUpper[iColumn]) {
983                    infeasibility += integerTolerance+1.0e-12;
984                    preferredWay=1;
985                  } else {
986                    infeasibility += integerTolerance+1.0e-12;
987                    preferredWay=1;
988                  }
989                } else {
990                  priorityLevel += 10000000;
991                }
992              } else {
993                // switch off if not possible
994                bestSolution=NULL;
995                model->setHotstartStrategy(0);
996              }
997            }
998          }
999        }
1000        if (infeasibility) {
1001          // Increase estimated degradation to solution
1002          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
1003          numberUnsatisfied_++;
1004          // Better priority? Flush choices.
1005          if (priorityLevel<bestPriority) {
1006            int j;
1007            iSmallest=0;
1008            for (j=0;j<maximumStrong;j++) {
1009              choice[j].upMovement=0.0;
1010              delete choice[j].possibleBranch;
1011              choice[j].possibleBranch=NULL;
1012            }
1013            bestPriority = priorityLevel;
1014            mostAway=1.0e-100;
1015            numberStrong=0;
1016          } else if (priorityLevel>bestPriority) {
1017            continue;
1018          }
1019          // Check for suitability based on infeasibility.
1020          if (infeasibility>mostAway) {
1021            //add to list
1022            choice[iSmallest].upMovement=infeasibility;
1023            delete choice[iSmallest].possibleBranch;
1024            choice[iSmallest].possibleBranch=object->createBranch(preferredWay);
1025            numberStrong = CoinMax(numberStrong,iSmallest+1);
1026            // Save which object it was
1027            choice[iSmallest].objectNumber=i;
1028            int j;
1029            iSmallest=-1;
1030            mostAway = 1.0e50;
1031            for (j=0;j<maximumStrong;j++) {
1032              if (choice[j].upMovement<mostAway) {
1033                mostAway=choice[j].upMovement;
1034                iSmallest=j;
1035              }
1036            }
1037          }
1038        }
1039      }
1040      if (numberUnsatisfied_) {
1041        // some infeasibilities - go to next steps
1042        break;
1043      } else if (!iPass) {
1044        // looks like a solution - get paranoid
1045        bool roundAgain=false;
1046        // get basis
1047        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1048        if (!ws)
1049          break;
1050        for (i=0;i<numberColumns;i++) {
1051          double value = saveSolution[i];
1052          if (value<lower[i]) {
1053            saveSolution[i]=lower[i];
1054            roundAgain=true;
1055            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1056          } else if (value>upper[i]) {
1057            saveSolution[i]=upper[i];
1058            roundAgain=true;
1059            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1060          } 
1061        }
1062        if (roundAgain) {
1063          // restore basis
1064          solver->setWarmStart(ws);
1065          delete ws;
1066          solver->resolve();
1067          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1068          model->reserveCurrentSolution(saveSolution);
1069          if (!solver->isProvenOptimal()) {
1070            // infeasible
1071            anyAction=-2;
1072            break;
1073          }
1074        } else {
1075          delete ws;
1076          break;
1077        }
1078      }
1079    }
1080    /* Some solvers can do the strong branching calculations faster if
1081       they do them all at once.  At present only Clp does for ordinary
1082       integers but I think this coding would be easy to modify
1083    */
1084    bool allNormal=true; // to say if we can do fast strong branching
1085    // Say which one will be best
1086    int bestChoice=0;
1087    double worstInfeasibility=0.0;
1088    for (i=0;i<numberStrong;i++) {
1089      choice[i].numIntInfeasUp = numberUnsatisfied_;
1090      choice[i].numIntInfeasDown = numberUnsatisfied_;
1091      choice[i].fix=0; // say not fixed
1092      if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1093        allNormal=false; // Something odd so lets skip clever fast branching
1094      if ( !model->object(choice[i].objectNumber)->boundBranch())
1095        numberStrong=0; // switch off
1096      // Do best choice in case switched off
1097      if (choice[i].upMovement>worstInfeasibility) {
1098        worstInfeasibility=choice[i].upMovement;
1099        bestChoice=i;
1100      }
1101    }
1102    // If we have hit max time don't do strong branching
1103    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1104                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1105    // also give up if we are looping round too much
1106    if (hitMaxTime||numberPassesLeft<=0)
1107      numberStrong=0;
1108    /*
1109      Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1110      fall through to simple branching.
1111     
1112      Setup for strong branching involves saving the current basis (for restoration
1113      afterwards) and setting up for hot starts.
1114    */
1115    if (numberStrong&&saveNumberStrong) {
1116     
1117      bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1118      solveAll=true;
1119      // worth trying if too many times
1120      // Save basis
1121      CoinWarmStart * ws = solver->getWarmStart();
1122      // save limit
1123      int saveLimit;
1124      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1125      if (beforeSolution&&saveLimit<100)
1126        solver->setIntParam(OsiMaxNumIterationHotStart,100); // go to end
1127     
1128      /* If we are doing all strong branching in one go then we create new arrays
1129         to store information.  If clp NULL then doing old way.
1130         Going down -
1131         outputSolution[2*i] is final solution.
1132         outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1133         outputStuff[2*i+numberStrong] is number iterations
1134         On entry newUpper[i] is new upper bound, on exit obj change
1135         Going up -
1136         outputSolution[2*i+1] is final solution.
1137         outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1138         outputStuff[2*i+1+numberStrong] is number iterations
1139       On entry newLower[i] is new lower bound, on exit obj change
1140      */
1141      OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1142      ClpSimplex * clp=NULL;
1143      double * newLower = NULL;
1144      double * newUpper = NULL;
1145      double ** outputSolution=NULL;
1146      int * outputStuff=NULL;
1147      // Go back to normal way if user wants it
1148      if (osiclp&&(osiclp->specialOptions()&16)!=0&&osiclp->specialOptions()>0)
1149        allNormal=false;
1150      if (osiclp&&!allNormal) {
1151        // say do fast
1152        int easy=1;
1153        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1154      }
1155      if (osiclp&& allNormal) {
1156        clp = osiclp->getModelPtr();
1157        // Clp - do a different way
1158        newLower = new double[numberStrong];
1159        newUpper = new double[numberStrong];
1160        outputSolution = new double * [2*numberStrong];
1161        outputStuff = new int [4*numberStrong];
1162        int * which = new int[numberStrong];
1163        int startFinishOptions;
1164        int specialOptions = osiclp->specialOptions();
1165        int clpOptions = clp->specialOptions();
1166        int returnCode=0;
1167#define CRUNCH
1168#ifdef CRUNCH
1169        // Crunch down problem
1170        int numberRows = clp->numberRows();
1171        // Use dual region
1172        double * rhs = clp->dualRowSolution();
1173        int * whichRow = new int[3*numberRows];
1174        int * whichColumn = new int[2*numberColumns];
1175        int nBound;
1176        ClpSimplex * small = ((ClpSimplexOther *) clp)->crunch(rhs,whichRow,whichColumn,nBound,true);
1177        if (!small) {
1178          anyAction=-2;
1179          //printf("XXXX Inf by inspection\n");
1180          delete [] whichColumn;
1181          whichColumn=NULL;
1182          delete [] whichRow;
1183          whichRow=NULL;
1184          break;
1185        } else {
1186          clp = small;
1187        }
1188#else
1189        int saveLogLevel = clp->logLevel();
1190        int saveMaxIts = clp->maximumIterations();
1191#endif
1192        clp->setLogLevel(0);
1193        if((specialOptions&1)==0) {
1194          startFinishOptions=0;
1195          clp->setSpecialOptions(clpOptions|(64|1024));
1196        } else {
1197          startFinishOptions=1+2+4;
1198          //startFinishOptions=1+4; // for moment re-factorize
1199          if((specialOptions&4)==0) 
1200            clp->setSpecialOptions(clpOptions|(64|128|512|1024|4096));
1201          else
1202            clp->setSpecialOptions(clpOptions|(64|128|512|1024|2048|4096));
1203        }
1204        // User may want to clean up before strong branching
1205        if ((clp->specialOptions()&32)!=0) {
1206          clp->primal(1);
1207          if (clp->numberIterations())
1208            model->messageHandler()->message(CBC_ITERATE_STRONG,model->messages())
1209              << clp->numberIterations()
1210              <<CoinMessageEol;
1211        }
1212        clp->setMaximumIterations(saveLimit);
1213#ifdef CRUNCH
1214        int * backColumn = whichColumn+numberColumns;
1215#endif
1216        for (i=0;i<numberStrong;i++) {
1217          int iObject = choice[i].objectNumber;
1218          const CbcObject * object = model->object(iObject);
1219          const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1220          int iSequence = simple->modelSequence();
1221          newLower[i]= ceil(saveSolution[iSequence]);
1222          newUpper[i]= floor(saveSolution[iSequence]);
1223#ifdef CRUNCH
1224          iSequence = backColumn[iSequence];
1225          assert (iSequence>=0);
1226#endif
1227          which[i]=iSequence;
1228          outputSolution[2*i]= new double [numberColumns];
1229          outputSolution[2*i+1]= new double [numberColumns];
1230        }
1231        //clp->writeMps("bad");
1232        returnCode=clp->strongBranching(numberStrong,which,
1233                                            newLower, newUpper,outputSolution,
1234                                            outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1235                                            startFinishOptions);
1236#ifndef CRUNCH
1237        clp->setSpecialOptions(clpOptions); // restore
1238        clp->setMaximumIterations(saveMaxIts);
1239        clp->setLogLevel(saveLogLevel);
1240#endif
1241        if (returnCode==-2) {
1242          // bad factorization!!!
1243          // Doing normal way
1244          // Mark hot start
1245          solver->markHotStart();
1246          clp = NULL;
1247        } else {
1248#ifdef CRUNCH
1249          // extract solution
1250          //bool checkSol=true;
1251          for (i=0;i<numberStrong;i++) {
1252            int iObject = choice[i].objectNumber;
1253            const CbcObject * object = model->object(iObject);
1254            const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1255            int iSequence = simple->modelSequence();
1256            which[i]=iSequence;
1257            double * sol = outputSolution[2*i];
1258            double * sol2 = outputSolution[2*i+1];
1259            //bool x=true;
1260            //bool x2=true;
1261            for (int iColumn=numberColumns-1;iColumn>=0;iColumn--) {
1262              int jColumn = backColumn[iColumn];
1263              if (jColumn>=0) {
1264                sol[iColumn]=sol[jColumn];
1265                sol2[iColumn]=sol2[jColumn];
1266              } else {
1267                sol[iColumn]=saveSolution[iColumn];
1268                sol2[iColumn]=saveSolution[iColumn];
1269              }
1270            }
1271          }
1272#endif
1273        }
1274#ifdef CRUNCH
1275        delete [] whichColumn;
1276        delete [] whichRow;
1277        delete small;
1278#endif
1279        delete [] which;
1280      } else {
1281        // Doing normal way
1282        // Mark hot start
1283        solver->markHotStart();
1284      }
1285      /*
1286        Open a loop to do the strong branching LPs. For each candidate variable,
1287        solve an LP with the variable forced down, then up. If a direction turns
1288        out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1289        force the objective change to be big (1.0e100). If we determine the problem
1290        is infeasible, or find a monotone variable, escape the loop.
1291       
1292        TODO: The `restore bounds' part might be better encapsulated as an
1293        unbranch() method. Branching objects more exotic than simple integers
1294        or cliques might not restrict themselves to variable bounds.
1295
1296        TODO: Virtuous solvers invalidate the current solution (or give bogus
1297        results :-) when the bounds are changed out from under them. So we
1298        need to do all the work associated with finding a new solution before
1299        restoring the bounds.
1300      */
1301      for (i = 0 ; i < numberStrong ; i++)
1302        { double objectiveChange ;
1303        double newObjectiveValue=1.0e100;
1304        // status is 0 finished, 1 infeasible and other
1305        int iStatus;
1306        /*
1307          Try the down direction first. (Specify the initial branching alternative as
1308          down with a call to way(-1). Each subsequent call to branch() performs the
1309          specified branch and advances the branch object state to the next branch
1310          alternative.)
1311        */
1312        if (!clp) {
1313          choice[i].possibleBranch->way(-1) ;
1314          choice[i].possibleBranch->branch() ;
1315          bool feasible=true;
1316          if (checkFeasibility) {
1317            // check branching did not make infeasible
1318            int iColumn;
1319            int numberColumns = solver->getNumCols();
1320            const double * columnLower = solver->getColLower();
1321            const double * columnUpper = solver->getColUpper();
1322            for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1323              if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1324                feasible=false;
1325            }
1326          }
1327          if (feasible) {
1328            solver->solveFromHotStart() ;
1329            /*
1330              We now have an estimate of objective degradation that we can use for strong
1331              branching. If we're over the cutoff, the variable is monotone up.
1332              If we actually made it to optimality, check for a solution, and if we have
1333              a good one, call setBestSolution to process it. Note that this may reduce the
1334              cutoff, so we check again to see if we can declare this variable monotone.
1335            */
1336            if (solver->isProvenOptimal())
1337              iStatus=0; // optimal
1338            else if (solver->isIterationLimitReached()
1339                     &&!solver->isDualObjectiveLimitReached())
1340              iStatus=2; // unknown
1341            else
1342              iStatus=1; // infeasible
1343            newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1344            choice[i].numItersDown = solver->getIterationCount();
1345          } else {
1346            iStatus=1; // infeasible
1347            newObjectiveValue = 1.0e100;
1348            choice[i].numItersDown = 0;
1349          }
1350          objectiveChange = newObjectiveValue-objectiveValue ;
1351        } else {
1352          iStatus = outputStuff[2*i];
1353          choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1354          newObjectiveValue = objectiveValue+newUpper[i];
1355          solver->setColSolution(outputSolution[2*i]);
1356        }
1357        objectiveChange = newObjectiveValue  - objectiveValue;
1358        if (!iStatus) {
1359          choice[i].finishedDown = true ;
1360          if (newObjectiveValue>=model->getCutoff()) {
1361            objectiveChange = 1.0e100; // say infeasible
1362          } else {
1363            // See if integer solution
1364            if (model->feasibleSolution(choice[i].numIntInfeasDown,
1365                                        choice[i].numObjInfeasDown)
1366                &&model->problemFeasibility()->feasible(model,-1)>=0) {
1367              model->setBestSolution(CBC_STRONGSOL,
1368                                     newObjectiveValue,
1369                                     solver->getColSolution()) ;
1370              model->setLastHeuristic(NULL);
1371              model->incrementUsed(solver->getColSolution());
1372              if (newObjectiveValue >= model->getCutoff())      //  *new* cutoff
1373                objectiveChange = 1.0e100 ;
1374            }
1375          }
1376        } else if (iStatus==1) {
1377          objectiveChange = 1.0e100 ;
1378        } else {
1379          // Can't say much as we did not finish
1380          choice[i].finishedDown = false ;
1381        }
1382        choice[i].downMovement = objectiveChange ;
1383       
1384        // restore bounds
1385        if (!clp)
1386          { for (int j=0;j<numberColumns;j++) {
1387            if (saveLower[j] != lower[j])
1388              solver->setColLower(j,saveLower[j]);
1389            if (saveUpper[j] != upper[j])
1390              solver->setColUpper(j,saveUpper[j]);
1391          }
1392          }
1393        //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1394        //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
1395        //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
1396        //     choice[i].numObjInfeasDown);
1397       
1398        // repeat the whole exercise, forcing the variable up
1399        if (!clp) {
1400          choice[i].possibleBranch->branch();
1401          bool feasible=true;
1402          if (checkFeasibility) {
1403            // check branching did not make infeasible
1404            int iColumn;
1405            int numberColumns = solver->getNumCols();
1406            const double * columnLower = solver->getColLower();
1407            const double * columnUpper = solver->getColUpper();
1408            for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1409              if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1410                feasible=false;
1411            }
1412          }
1413          if (feasible) {
1414            solver->solveFromHotStart() ;
1415            /*
1416              We now have an estimate of objective degradation that we can use for strong
1417              branching. If we're over the cutoff, the variable is monotone up.
1418              If we actually made it to optimality, check for a solution, and if we have
1419              a good one, call setBestSolution to process it. Note that this may reduce the
1420              cutoff, so we check again to see if we can declare this variable monotone.
1421            */
1422            if (solver->isProvenOptimal())
1423              iStatus=0; // optimal
1424            else if (solver->isIterationLimitReached()
1425                     &&!solver->isDualObjectiveLimitReached())
1426              iStatus=2; // unknown
1427            else
1428              iStatus=1; // infeasible
1429            newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1430            choice[i].numItersUp = solver->getIterationCount();
1431          } else {
1432            iStatus=1; // infeasible
1433            newObjectiveValue = 1.0e100;
1434            choice[i].numItersDown = 0;
1435          }
1436          objectiveChange = newObjectiveValue-objectiveValue ;
1437        } else {
1438          iStatus = outputStuff[2*i+1];
1439          choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
1440          newObjectiveValue = objectiveValue+newLower[i];
1441          solver->setColSolution(outputSolution[2*i+1]);
1442        }
1443        objectiveChange = newObjectiveValue  - objectiveValue;
1444        if (!iStatus) {
1445          choice[i].finishedUp = true ;
1446          if (newObjectiveValue>=model->getCutoff()) {
1447            objectiveChange = 1.0e100; // say infeasible
1448          } else {
1449            // See if integer solution
1450            if (model->feasibleSolution(choice[i].numIntInfeasUp,
1451                                        choice[i].numObjInfeasUp)
1452                &&model->problemFeasibility()->feasible(model,-1)>=0) {
1453              model->setBestSolution(CBC_STRONGSOL,
1454                                     newObjectiveValue,
1455                                     solver->getColSolution()) ;
1456              model->setLastHeuristic(NULL);
1457              model->incrementUsed(solver->getColSolution());
1458              if (newObjectiveValue >= model->getCutoff())      //  *new* cutoff
1459                objectiveChange = 1.0e100 ;
1460            }
1461          }
1462        } else if (iStatus==1) {
1463          objectiveChange = 1.0e100 ;
1464        } else {
1465          // Can't say much as we did not finish
1466          choice[i].finishedUp = false ;
1467        }
1468        choice[i].upMovement = objectiveChange ;
1469       
1470        // restore bounds
1471        if (!clp)
1472          { for (int j=0;j<numberColumns;j++) {
1473            if (saveLower[j] != lower[j])
1474              solver->setColLower(j,saveLower[j]);
1475            if (saveUpper[j] != upper[j])
1476              solver->setColUpper(j,saveUpper[j]);
1477          }
1478          }
1479       
1480        //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1481        //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
1482        //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
1483        //     choice[i].numObjInfeasUp);
1484       
1485        /*
1486          End of evaluation for this candidate variable. Possibilities are:
1487          * Both sides below cutoff; this variable is a candidate for branching.
1488          * Both sides infeasible or above the objective cutoff: no further action
1489          here. Break from the evaluation loop and assume the node will be purged
1490          by the caller.
1491          * One side below cutoff: Install the branch (i.e., fix the variable). Break
1492          from the evaluation loop and assume the node will be reoptimised by the
1493          caller.
1494        */
1495        if (choice[i].upMovement<1.0e100) {
1496          if(choice[i].downMovement<1.0e100) {
1497            // feasible - no action
1498          } else {
1499            // up feasible, down infeasible
1500            anyAction=-1;
1501            //printf("Down infeasible for choice %d sequence %d\n",i,
1502            // model->object(choice[i].objectNumber)->columnNumber());
1503            if (!solveAll) {
1504              choice[i].possibleBranch->way(1);
1505              choice[i].possibleBranch->branch();
1506              break;
1507            } else {
1508              choice[i].fix=1;
1509            }
1510          }
1511        } else {
1512          if(choice[i].downMovement<1.0e100) {
1513            // down feasible, up infeasible
1514            anyAction=-1;
1515            //printf("Up infeasible for choice %d sequence %d\n",i,
1516            // model->object(choice[i].objectNumber)->columnNumber());
1517            if (!solveAll) {
1518              choice[i].possibleBranch->way(-1);
1519              choice[i].possibleBranch->branch();
1520              break;
1521            } else {
1522              choice[i].fix=-1;
1523            }
1524          } else {
1525            // neither side feasible
1526            anyAction=-2;
1527            //printf("Both infeasible for choice %d sequence %d\n",i,
1528            // model->object(choice[i].objectNumber)->columnNumber());
1529            break;
1530          }
1531        }
1532        bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1533                            model->getDblParam(CbcModel::CbcMaximumSeconds));
1534        if (hitMaxTime) {
1535          numberStrong=i+1;
1536          break;
1537        }
1538        }
1539      if (!clp) {
1540        // Delete the snapshot
1541        solver->unmarkHotStart();
1542      } else {
1543        delete [] newLower;
1544        delete [] newUpper;
1545        delete [] outputStuff;
1546        int i;
1547        for (i=0;i<2*numberStrong;i++)
1548          delete [] outputSolution[i];
1549        delete [] outputSolution;
1550      }
1551      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
1552      // restore basis
1553      solver->setWarmStart(ws);
1554      // Unless infeasible we will carry on
1555      // But we could fix anyway
1556      if (anyAction==-1&&solveAll) {
1557        // apply and take off
1558        for (i = 0 ; i < numberStrong ; i++) {
1559          if (choice[i].fix) {
1560            choice[i].possibleBranch->way(choice[i].fix) ;
1561            choice[i].possibleBranch->branch() ;
1562          }
1563        }
1564        bool feasible=true;
1565        if (checkFeasibility) {
1566          // check branching did not make infeasible
1567          int iColumn;
1568          int numberColumns = solver->getNumCols();
1569          const double * columnLower = solver->getColLower();
1570          const double * columnUpper = solver->getColUpper();
1571          for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1572            if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1573              feasible=false;
1574          }
1575        }
1576        if (feasible) {
1577          // can do quick optimality check
1578          int easy=2;
1579          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1580          solver->resolve() ;
1581          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1582          feasible = solver->isProvenOptimal();
1583        }
1584        if (feasible) {
1585          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1586          model->reserveCurrentSolution(saveSolution);
1587          memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
1588          memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
1589          // Clean up all candidates whih are fixed
1590          int numberLeft=0;
1591          for (i = 0 ; i < numberStrong ; i++) {
1592            CbcStrongInfo thisChoice = choice[i];
1593            choice[i].possibleBranch=NULL;
1594            const CbcObject * object = model->object(thisChoice.objectNumber);
1595            int preferredWay;
1596            double infeasibility = object->infeasibility(preferredWay);
1597            if (!infeasibility) {
1598              // take out
1599              delete thisChoice.possibleBranch;
1600            } else {
1601              choice[numberLeft++]=thisChoice;
1602            }
1603          }
1604          numberStrong=numberLeft;
1605          for (;i<maximumStrong;i++) {
1606            delete choice[i].possibleBranch;
1607            choice[i].possibleBranch=NULL;
1608          }
1609          // If all fixed then round again
1610          if (!numberLeft) {
1611            finished=false;
1612            numberStrong=0;
1613            saveNumberStrong=0;
1614            maximumStrong=1;
1615          } else {
1616            anyAction=0;
1617          }
1618          // If these two uncommented then different action
1619          anyAction=-1;
1620          finished=true;
1621          //printf("some fixed but continuing %d left\n",numberLeft);
1622        } else {
1623          anyAction=-2; // say infeasible
1624        }
1625      }
1626      delete ws;
1627     
1628      /*
1629        anyAction >= 0 indicates that strong branching didn't produce any monotone
1630        variables. Sift through the candidates for the best one.
1631       
1632        QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
1633        local to this code block. Perhaps should be numberNodes_ from model?
1634        Unclear what this calculation is doing.
1635      */
1636      if (anyAction>=0) {
1637       
1638        int numberNodes = model->getNodeCount();
1639        // get average cost per iteration and assume stopped ones
1640        // would stop after 50% more iterations at average cost??? !!! ???
1641        double averageCostPerIteration=0.0;
1642        double totalNumberIterations=1.0;
1643        int smallestNumberInfeasibilities=INT_MAX;
1644        for (i=0;i<numberStrong;i++) {
1645          totalNumberIterations += choice[i].numItersDown +
1646            choice[i].numItersUp ;
1647          averageCostPerIteration += choice[i].downMovement +
1648            choice[i].upMovement;
1649          smallestNumberInfeasibilities= 
1650            CoinMin(CoinMin(choice[i].numIntInfeasDown ,
1651                            choice[i].numIntInfeasUp ),
1652                    smallestNumberInfeasibilities);
1653        }
1654        //if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
1655        //numberNodes=1000000; // switch off search for better solution
1656        numberNodes=1000000; // switch off anyway
1657        averageCostPerIteration /= totalNumberIterations;
1658        // all feasible - choose best bet
1659       
1660        // New method does all at once so it can be more sophisticated
1661        // in deciding how to balance actions.
1662        // But it does need arrays
1663        double * changeUp = new double [numberStrong];
1664        int * numberInfeasibilitiesUp = new int [numberStrong];
1665        double * changeDown = new double [numberStrong];
1666        int * numberInfeasibilitiesDown = new int [numberStrong];
1667        CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
1668        for (i = 0 ; i < numberStrong ; i++) {
1669          int iColumn = choice[i].possibleBranch->variable() ;
1670          model->messageHandler()->message(CBC_STRONG,model->messages())
1671            << i << iColumn
1672            <<choice[i].downMovement<<choice[i].numIntInfeasDown
1673            <<choice[i].upMovement<<choice[i].numIntInfeasUp
1674            <<choice[i].possibleBranch->value()
1675            <<CoinMessageEol;
1676          changeUp[i]=choice[i].upMovement;
1677          numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1678          changeDown[i]=choice[i].downMovement;
1679          numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1680          objects[i] = choice[i].possibleBranch;
1681        }
1682        int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
1683                                               changeUp,numberInfeasibilitiesUp,
1684                                               changeDown,numberInfeasibilitiesDown,
1685                                               objectiveValue);
1686        // move branching object and make sure it will not be deleted
1687        if (whichObject>=0) {
1688          branch_ = objects[whichObject];
1689          choice[whichObject].possibleBranch=NULL;
1690        }
1691        delete [] changeUp;
1692        delete [] numberInfeasibilitiesUp;
1693        delete [] changeDown;
1694        delete [] numberInfeasibilitiesDown;
1695        delete [] objects;
1696      }
1697      if (osiclp&&!allNormal) {
1698        // back to normal
1699        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1700      }
1701    }
1702    /*
1703      Simple branching. Probably just one, but we may have got here
1704      because of an odd branch e.g. a cut
1705    */
1706    else {
1707      // not strong
1708      // C) create branching object
1709      branch_ = choice[bestChoice].possibleBranch;
1710      choice[bestChoice].possibleBranch=NULL;
1711    }
1712  }
1713  // Set guessed solution value
1714  objectiveValue_ = solver->getObjSense()*saveObjectiveValue;
1715  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
1716/*
1717  Cleanup, then we're outta here.
1718*/
1719  if (!model->branchingMethod())
1720    delete decision;
1721   
1722  for (i=0;i<maximumStrong;i++)
1723    delete choice[i].possibleBranch;
1724  delete [] choice;
1725  delete [] saveLower;
1726  delete [] saveUpper;
1727 
1728  // restore solution
1729  solver->setColSolution(saveSolution);
1730  delete [] saveSolution;
1731  return anyAction;
1732}
1733
1734/*
1735  Version for dynamic pseudo costs.
1736
1737  **** For now just return if anything odd
1738  later allow even if odd
1739
1740  The routine scans through the object list of the model looking for objects
1741  that indicate infeasibility. It tests each object using strong branching
1742  and selects the one with the least objective degradation.  A corresponding
1743  branching object is left attached to lastNode.
1744  This version gives preference in evaluation to variables which
1745  have not been evaluated many times.  It also uses numberStrong
1746  to say give up if last few tries have not changed incumbent.
1747  See Achterberg, Koch and Martin.
1748
1749  If strong branching is disabled, a candidate object is chosen essentially
1750  at random (whatever object ends up in pos'n 0 of the candidate array).
1751
1752  If a branching candidate is found to be monotone, bounds are set to fix the
1753  variable and the routine immediately returns (the caller is expected to
1754  reoptimize).
1755
1756  If a branching candidate is found to result in infeasibility in both
1757  directions, the routine immediately returns an indication of infeasibility.
1758
1759  Returns:  0   both branch directions are feasible
1760           -1   branching variable is monotone
1761           -2   infeasible
1762           -3   Use another method
1763
1764           For now just fix on objective from strong branching.
1765*/
1766
1767int CbcNode::chooseDynamicBranch (CbcModel *model, CbcNode *lastNode,int numberPassesLeft)
1768
1769{ if (lastNode)
1770    depth_ = lastNode->depth_+1;
1771  else
1772    depth_ = 0;
1773  delete branch_;
1774  branch_=NULL;
1775  OsiSolverInterface * solver = model->solver();
1776  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
1777  double cutoff =model->getCutoff();
1778  double distanceToCutoff = cutoff-objectiveValue_;
1779  const double * lower = solver->getColLower();
1780  const double * upper = solver->getColUpper();
1781  // See what user thinks
1782  int anyAction=model->problemFeasibility()->feasible(model,0);
1783  if (anyAction) {
1784    // will return -2 if infeasible , 0 if treat as integer
1785    return anyAction-1;
1786  }
1787  int i;
1788  int stateOfSearch = model->stateOfSearch();
1789  int numberStrong=model->numberStrong();
1790  // But make more likely to get out after some times
1791  int changeStrategy=numberStrong;
1792  double changeFactor=1.0;
1793  // Use minimum of this and one stored in objects
1794  //int numberBeforeTrust = model->numberBeforeTrust();
1795  int numberObjects = model->numberObjects();
1796  bool checkFeasibility = numberObjects>model->numberIntegers();
1797  // For now return if not simple
1798  if (checkFeasibility)
1799    return -3;
1800  // Return if doing hot start (in BAB sense)
1801  int hotstartStrategy=model->getHotstartStrategy();
1802  if (hotstartStrategy>0) 
1803    return -3;
1804  // Pass number
1805  int kPass=0;
1806  int numberColumns = model->getNumCols();
1807  int numberRows = solver->getNumRows();
1808  double * saveUpper = new double[numberColumns];
1809  double * saveLower = new double[numberColumns];
1810  // Ratio to cutoff for pseudo costs
1811  double tryStrongPseudo = 0.8*distanceToCutoff;
1812  // Ratio to cutoff for penalties
1813  //double tryStrongPenalty = 0.5*distanceToCutoff;
1814
1815  // Save solution in case heuristics need good solution later
1816 
1817  double * saveSolution = new double[numberColumns];
1818  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1819  model->reserveCurrentSolution(saveSolution);
1820  /*
1821    Get a branching decision object. Use the default dynamic decision criteria unless
1822    the user has loaded a decision method into the model.
1823  */
1824  CbcBranchDecision *decision = model->branchingMethod();
1825  if (!decision)
1826    decision = new CbcBranchDynamicDecision();
1827
1828  for (i=0;i<numberColumns;i++) {
1829    saveLower[i] = lower[i];
1830    saveUpper[i] = upper[i];
1831  }
1832  // Get arrays to sort
1833  double * sort = new double[numberObjects];
1834  int * whichObject = new int[numberObjects];
1835  int * objectMark = new int[2*numberObjects+1];
1836  // Arrays with movements
1837  double * upEstimate = new double[numberObjects];
1838  double * downEstimate = new double[numberObjects];
1839  CbcStrongInfo * fixObject = new CbcStrongInfo[numberObjects];
1840  double estimatedDegradation=0.0; 
1841  int numberNodes=model->getNodeCount();
1842  int numberBeforeTrust = model->numberBeforeTrust();
1843  int numberPenalties = model->numberPenalties();
1844  if (numberBeforeTrust>=1000000) {
1845    numberBeforeTrust = numberBeforeTrust % 1000000;
1846    numberPenalties=0;
1847  } else if (numberBeforeTrust<0) {
1848    numberPenalties=numberColumns;
1849    numberBeforeTrust=0;
1850  }
1851  double scaleFactor = model->penaltyScaleFactor();
1852  // May go round twice if strong branching fixes all local candidates
1853  bool finished=false;
1854  while(!finished) {
1855    finished=true;
1856    decision->initialize(model);
1857    // Some objects may compute an estimate of best solution from here
1858    estimatedDegradation=0.0; 
1859    int numberToFix=0;
1860    int numberIntegerInfeasibilities=0; // without odd ones
1861    int numberToDo=0;
1862    double averageDown=0.0;
1863    int numberDown=0;
1864    double averageUp=0.0;
1865    int numberUp=0;
1866    int iBestNot=-1;
1867    int iBestGot=-1;
1868    double best=0.0;
1869    int numberNotTrusted=0;
1870   
1871    // We may go round this loop twice (only if we think we have solution)
1872    for (int iPass=0;iPass<2;iPass++) {
1873     
1874      // compute current state
1875      int numberObjectInfeasibilities; // just odd ones
1876      model->feasibleSolution(
1877                              numberIntegerInfeasibilities,
1878                              numberObjectInfeasibilities);
1879     
1880      // Some objects may compute an estimate of best solution from here
1881      estimatedDegradation=0.0; 
1882      numberUnsatisfied_ = 0;
1883      int bestPriority=INT_MAX;
1884      /*
1885        Scan for branching objects that indicate infeasibility. Choose candidates
1886        using priority as the first criteria, then integer infeasibility.
1887       
1888        The algorithm is to fill the array with a set of good candidates (by
1889        infeasibility) with priority bestPriority.  Finding a candidate with
1890        priority better (less) than bestPriority flushes the choice array. (This
1891        serves as initialization when the first candidate is found.)
1892       
1893      */
1894      numberToDo=0;
1895      averageDown=0.0;
1896      numberDown=0;
1897      averageUp=0.0;
1898      numberUp=0;
1899      iBestNot=-1;
1900      double bestNot=0.0;
1901      iBestGot=-1;
1902      best=0.0;
1903#define PRINT_STUFF -1
1904      for (i=0;i<numberObjects;i++) {
1905        CbcObject * object = model->modifiableObject(i);
1906        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
1907          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
1908        assert(dynamicObject);
1909        int preferredWay;
1910        double infeasibility = object->infeasibility(preferredWay);
1911        int priorityLevel = object->priority();
1912        bool gotDown=false;
1913        int numberThisDown = dynamicObject->numberTimesDown();
1914        if (numberThisDown) {
1915          averageDown += dynamicObject->downDynamicPseudoCost();
1916          numberDown++;
1917          if (numberThisDown>=numberBeforeTrust)
1918            gotDown=true;
1919        }
1920        bool gotUp=false;
1921        int numberThisUp = dynamicObject->numberTimesUp();
1922        if (numberThisUp) {
1923          averageUp += dynamicObject->upDynamicPseudoCost();
1924          numberUp++;
1925          if (numberThisUp>=numberBeforeTrust)
1926            gotUp=true;
1927        }
1928        if (infeasibility) {
1929          if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
1930            printf("%d down %d %g up %d %g - infeas %g\n",
1931                   i,numberThisDown,object->downEstimate(),numberThisUp,object->upEstimate(),
1932                   infeasibility);
1933          // Increase estimated degradation to solution
1934          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
1935          downEstimate[i]=object->downEstimate();
1936          upEstimate[i]=object->upEstimate();
1937          numberUnsatisfied_++;
1938          // Better priority? Flush choices.
1939          if (priorityLevel<bestPriority) {
1940            numberToDo=0;
1941            bestPriority = priorityLevel;
1942            iBestGot=-1;
1943            best=0.0;
1944            numberNotTrusted=0;
1945          } else if (priorityLevel>bestPriority) {
1946            continue;
1947          }
1948          if (!gotUp||!gotDown)
1949            numberNotTrusted++;
1950          // Check for suitability based on infeasibility.
1951          if ((gotDown&&gotUp)&&numberStrong>0) {
1952            sort[numberToDo]=-infeasibility;
1953            if (infeasibility>best) {
1954              best=infeasibility;
1955              iBestGot=numberToDo;
1956            }
1957            if (infeasibility<tryStrongPseudo)
1958              objectMark[i]=2*numberBeforeTrust; // range info not wanted
1959            else
1960              objectMark[i]=numberBeforeTrust; // keep as possible cutoff
1961          } else {
1962            int iColumn = dynamicObject->columnNumber();
1963            double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
1964            sort[numberToDo]=part;
1965            if (1.0-fabs(part-0.5)>bestNot) {
1966              iBestNot=numberToDo;
1967              bestNot = 1.0-fabs(part-0.5);
1968            }
1969            objectMark[i]=numberThisDown+numberThisUp;
1970          }
1971          whichObject[numberToDo++]=i;
1972        } else {
1973          // for debug
1974          downEstimate[i]=-1.0;
1975          upEstimate[i]=-1.0;
1976        }
1977      }
1978      if (numberUnsatisfied_) {
1979        // some infeasibilities - go to next steps
1980        break;
1981      } else if (!iPass) {
1982        // looks like a solution - get paranoid
1983        bool roundAgain=false;
1984        // get basis
1985        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1986        if (!ws)
1987          break;
1988        double tolerance;
1989        solver->getDblParam(OsiPrimalTolerance,tolerance);
1990        for (i=0;i<numberColumns;i++) {
1991          double value = saveSolution[i];
1992          if (value<lower[i]-tolerance) {
1993            saveSolution[i]=lower[i];
1994            roundAgain=true;
1995            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1996          } else if (value>upper[i]+tolerance) {
1997            saveSolution[i]=upper[i];
1998            roundAgain=true;
1999            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
2000          } 
2001        }
2002        if (roundAgain) {
2003          // restore basis
2004          solver->setWarmStart(ws);
2005          solver->setColSolution(saveSolution);
2006          delete ws;
2007          bool takeHint;
2008          OsiHintStrength strength;
2009          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
2010          solver->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
2011          solver->resolve();
2012          solver->setHintParam(OsiDoDualInResolve,takeHint,strength) ;
2013          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2014          model->reserveCurrentSolution(saveSolution);
2015          if (!solver->isProvenOptimal()) {
2016            // infeasible
2017            anyAction=-2;
2018            break;
2019          }
2020        } else {
2021          delete ws;
2022          break;
2023        }
2024      }
2025    }
2026    if (anyAction==-2) {
2027      break;
2028    }
2029    bool solveAll=false; // set true to say look at all even if some fixed (experiment)
2030    solveAll=true;
2031    // skip if solution
2032    if (!numberUnsatisfied_)
2033      break;
2034    bool skipAll = (numberBeforeTrust>20&&numberNodes>20000&&numberNotTrusted==0);
2035    if ((numberNodes%20)==0)
2036      skipAll=false;
2037    // worth trying if too many times
2038    // Save basis
2039    CoinWarmStart * ws = NULL;
2040    // save limit
2041    int saveLimit=0;
2042    if (!skipAll) {
2043      ws = solver->getWarmStart();
2044      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
2045      if (!stateOfSearch&&saveLimit<100)
2046        solver->setIntParam(OsiMaxNumIterationHotStart,100); 
2047    }
2048    // Say which one will be best
2049    int whichChoice=0;
2050    int bestChoice;
2051    if (iBestGot>=0)
2052      bestChoice=iBestGot;
2053    else
2054      bestChoice=iBestNot;
2055    assert (bestChoice>=0);
2056    // If we have hit max time don't do strong branching
2057    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2058                        model->getDblParam(CbcModel::CbcMaximumSeconds));
2059    // also give up if we are looping round too much
2060    if (hitMaxTime||numberPassesLeft<=0||(!numberNotTrusted&&false)) {
2061      int iObject = whichObject[bestChoice];
2062      CbcObject * object = model->modifiableObject(iObject);
2063      int preferredWay;
2064      object->infeasibility(preferredWay);
2065      branch_=object->createBranch(preferredWay);
2066      branch_->way(preferredWay);
2067      delete ws;
2068      ws=NULL;
2069      break;
2070    } else {
2071      // say do fast
2072      int easy=1;
2073      if (!skipAll)
2074        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
2075      if (numberDown)
2076        averageDown /= (double) numberDown;
2077      else
2078        averageDown=1.0;
2079      if (numberUp)
2080        averageUp /= (double) numberUp;
2081      else
2082        averageUp=1.0;
2083      double average = 0.5*(averageUp+averageDown);
2084      int iDo;
2085      // Bias by trust
2086      int iObj;
2087      for (iObj=0;iObj<numberToDo;iObj++) {
2088        int iObject=whichObject[iObj];
2089        double add = objectMark[iObject];
2090        if (sort[iObj]<0.0)
2091          sort[iObj] += add*average;
2092        else
2093          sort[iObj] = -(sort[iObj]*averageDown+(1.0-sort[iObj])*averageUp);
2094      }
2095      // Sort
2096      CoinSort_2(sort,sort+numberToDo,whichObject);
2097      int needed2=numberToDo;
2098#define RANGING
2099#ifdef RANGING
2100      OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
2101      if (skipAll&&numberBeforeTrust)
2102        numberPenalties=0;
2103      if (osiclp&&numberPenalties) {
2104        // just get those not touched and best and not trusted
2105        int n=CoinMin(numberPenalties,numberToDo);
2106        int * which = objectMark+numberObjects+1;
2107        int i;
2108        int needed=0;
2109        for ( i=0;i<n;i++) {
2110          int iObject=whichObject[i];
2111          CbcObject * object = model->modifiableObject(iObject);
2112          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2113            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2114          which[needed++]=dynamicObject->columnNumber();
2115        }
2116        which--;
2117        which[0]=needed;
2118        assert (needed);
2119        osiclp->passInRanges(which);
2120        // Mark hot start and get ranges
2121        if (kPass) {
2122          // until can work out why solution can go funny
2123          int save = osiclp->specialOptions();
2124          osiclp->setSpecialOptions(save|256);
2125          solver->markHotStart();
2126          osiclp->setSpecialOptions(save);
2127        } else {
2128          solver->markHotStart();
2129        }
2130        kPass++;
2131        osiclp->passInRanges(NULL);
2132        needed2=0;
2133        int put=0;
2134        const double * downCost=osiclp->upRange();
2135        const double * upCost=osiclp->downRange();
2136        // Bug - so switch off for now
2137        double distanceToCutoff=COIN_DBL_MAX;
2138        //printf("numberTodo %d needed %d numberpenalties %d\n",numberToDo,needed,numberPenalties);
2139        for ( i=0;i<numberToDo;i++) {
2140          int iObject = whichObject[i];
2141          CbcObject * object = model->modifiableObject(iObject);
2142          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2143            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2144          int iSequence=dynamicObject->columnNumber();
2145          double estimate = sort[i];
2146          if (i<needed) {
2147            // We have looked at penalties
2148            int iAction=0;
2149            double value = saveSolution[iSequence];
2150            value -= floor(value);
2151            double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
2152            double downPenalty = CoinMin(downCost[i],1.0e110)*value;
2153            if (!numberBeforeTrust) {
2154              // override
2155              downEstimate[iObject]=downPenalty;
2156              upEstimate[iObject]=upPenalty;
2157            }
2158            if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
2159              printf("%d pen down ps %g -> %g up ps %g -> %g\n",
2160                     iObject,downCost[i],downPenalty,upCost[i],upPenalty);
2161            if (upPenalty>distanceToCutoff) {
2162              if(downPenalty>distanceToCutoff) {
2163                //printf("%d infeas both penalty %g %g\n",iObject,upPenalty,downPenalty);
2164                iAction=3;
2165              } else {
2166                //printf("%d infeas up penalty %g %g\n",iObject,upPenalty,downPenalty);
2167                iAction=2;
2168              }
2169            } else if(downPenalty>distanceToCutoff) {
2170              //printf("%d infeas down penalty %g %g\n",iObject,upPenalty,downPenalty);
2171              iAction=1;
2172            }
2173            if (iAction) {
2174              if (iAction==3) {
2175                //printf("%d infeas both penalty %g %g\n",iObject,upPenalty,downPenalty);
2176                anyAction=-2;
2177                delete ws;
2178                ws=NULL;
2179                break;
2180              } else if (iAction==2) {
2181                //printf("%d infeas up penalty %g %g\n",iObject,upPenalty,downPenalty);
2182                CbcStrongInfo choice;
2183                choice.objectNumber=iObject;
2184                choice.fix=-1;
2185                choice.possibleBranch=object->createBranch(-1);
2186                fixObject[numberToFix++]=choice;
2187                if (!anyAction)
2188                  anyAction=-1;
2189              } else {
2190                //printf("%d infeas down penalty %g %g\n",iObject,upPenalty,downPenalty);
2191                CbcStrongInfo choice;
2192                choice.objectNumber=iObject;
2193                choice.fix=1;
2194                choice.possibleBranch=object->createBranch(1);
2195                fixObject[numberToFix++]=choice;
2196                if (!anyAction)
2197                  anyAction=-1;
2198              }
2199            } else {
2200              double add = objectMark[iObject];
2201              double trueEstimate = sort[i] - add*average; 
2202              trueEstimate=0.0; // temp
2203              estimate = trueEstimate-scaleFactor*(upPenalty+downPenalty);
2204              sort[put]=estimate;
2205              whichObject[put++]=iObject;
2206              needed2=put;
2207            }
2208          } else {
2209            // just estimate
2210            double add = objectMark[iObject];
2211            double trueEstimate = sort[i] - add*average; 
2212            sort[put]=trueEstimate;
2213            whichObject[put++]=iObject;
2214          }
2215        }
2216        numberToDo=put;
2217      } else {
2218        if (!skipAll) {
2219          // Mark hot start
2220          solver->markHotStart();
2221          if (solver->isProvenPrimalInfeasible())
2222            printf("**** Hot start says node infeasible\n");
2223        }
2224        // make sure best will be first
2225        if (iBestGot>=0)
2226          sort[iBestGot]=-1.0e120;
2227      }
2228#else
2229      // Mark hot start
2230      solver->markHotStart();
2231      // make sure best will be first
2232      if (iBestGot>=0)
2233        sort[iBestGot]=-COIN_DBL_MAX;
2234#endif
2235      // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
2236#define ACTION 0
2237#if ACTION<2
2238      if (anyAction)
2239        numberToDo=0; // skip as we will be trying again
2240#endif
2241      // Sort (but just those re-computed)
2242      CoinSort_2(sort,sort+needed2,whichObject);
2243      // Just first if strong off
2244      if (!numberStrong)
2245        numberToDo=CoinMin(numberToDo,1);
2246      iDo=0;
2247      int saveLimit2;
2248      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit2);
2249      bool doQuickly = false; // numberToDo>2*numberStrong;
2250      //printf("todo %d, strong %d\n",numberToDo,numberStrong);
2251      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
2252      int numberTest2 = CoinMax(2*numberStrong,10)+1000000;
2253#if 0
2254      if (numberBeforeTrust>20&&(numberNodes>20000||(numberNodes>200&&numberNotTrusted==0))) {
2255        if ((numberNodes%20)!=0) {
2256          numberTest=0;
2257          doQuickly=true;
2258        }
2259      }
2260#else
2261      // Try nearly always off
2262      if ((numberNodes%20)!=0) {
2263        numberTest=0;
2264        doQuickly=true;
2265      } else {
2266        doQuickly=true;
2267        numberTest=2*numberStrong;
2268      }
2269#endif
2270      // if too many and big then just do 10 its
2271      if (!skipAll&&stateOfSearch) {
2272        if (numberNotTrusted>3*numberStrong&&numberRows>250&&numberColumns>1000)
2273          solver->setIntParam(OsiMaxNumIterationHotStart,10); 
2274      }
2275      double distanceToCutoff=model->getCutoff()-objectiveValue_;
2276      // larger and make negative for test
2277      distanceToCutoff *= -100.0;
2278      if (skipAll)
2279        distanceToCutoff = -COIN_DBL_MAX;
2280      for ( iDo=0;iDo<numberToDo;iDo++) {
2281        CbcStrongInfo choice;
2282        int iObject = whichObject[iDo];
2283        CbcObject * object = model->modifiableObject(iObject);
2284        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2285          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2286        int iColumn = dynamicObject->columnNumber();
2287        int preferredWay;
2288        object->infeasibility(preferredWay);
2289        choice.possibleBranch=object->createBranch(preferredWay);
2290        // Save which object it was
2291        choice.objectNumber=iObject;
2292        choice.numIntInfeasUp = numberUnsatisfied_;
2293        choice.numIntInfeasDown = numberUnsatisfied_;
2294        choice.upMovement = upEstimate[iObject];
2295        choice.downMovement = downEstimate[iObject];
2296        assert (choice.upMovement>=0.0);
2297        assert (choice.downMovement>=0.0);
2298        choice.fix=0; // say not fixed
2299        // see if can skip strong branching
2300        int canSkip = choice.possibleBranch->fillStrongInfo(choice);
2301        if (!doQuickly||numberTest>0)
2302          canSkip=0;
2303        if (!numberBeforeTrust) {
2304          choice.upMovement = upEstimate[iObject];
2305          choice.downMovement = downEstimate[iObject];
2306          canSkip=1;
2307        }
2308        if (sort[iDo]<distanceToCutoff)
2309          canSkip=0;
2310        if (numberTest2<=0&&sort[iDo]>distanceToCutoff)
2311          canSkip=1; // always skip
2312        if (model->messageHandler()->logLevel()>3&&numberBeforeTrust) 
2313          dynamicObject->print(1,choice.possibleBranch->value());
2314        // was if (!canSkip)
2315        if (!canSkip) {
2316          numberTest--;
2317          numberTest2--;
2318          // just do a few
2319          if (canSkip)
2320            solver->setIntParam(OsiMaxNumIterationHotStart,10); 
2321          double objectiveChange ;
2322          double newObjectiveValue=1.0e100;
2323          int j;
2324          // status is 0 finished, 1 infeasible and other
2325          int iStatus;
2326          /*
2327            Try the down direction first. (Specify the initial branching alternative as
2328            down with a call to way(-1). Each subsequent call to branch() performs the
2329            specified branch and advances the branch object state to the next branch
2330            alternative.)
2331          */
2332          choice.possibleBranch->way(-1) ;
2333          decision->saveBranchingObject( choice.possibleBranch);
2334          choice.possibleBranch->branch() ;
2335          solver->solveFromHotStart() ;
2336          /*
2337            We now have an estimate of objective degradation that we can use for strong
2338            branching. If we're over the cutoff, the variable is monotone up.
2339            If we actually made it to optimality, check for a solution, and if we have
2340            a good one, call setBestSolution to process it. Note that this may reduce the
2341            cutoff, so we check again to see if we can declare this variable monotone.
2342          */
2343          if (solver->isProvenOptimal())
2344            iStatus=0; // optimal
2345          else if (solver->isIterationLimitReached()
2346                   &&!solver->isDualObjectiveLimitReached())
2347            iStatus=2; // unknown
2348          else
2349            iStatus=1; // infeasible
2350          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2351          choice.numItersDown = solver->getIterationCount();
2352          objectiveChange = newObjectiveValue  - objectiveValue_;
2353          decision->updateInformation( solver,this);
2354          if (!iStatus) {
2355            choice.finishedDown = true ;
2356            if (newObjectiveValue>=cutoff) {
2357              objectiveChange = 1.0e100; // say infeasible
2358            } else {
2359              // See if integer solution
2360              if (model->feasibleSolution(choice.numIntInfeasDown,
2361                                          choice.numObjInfeasDown)
2362                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2363                model->setBestSolution(CBC_STRONGSOL,
2364                                       newObjectiveValue,
2365                                       solver->getColSolution()) ;
2366                model->setLastHeuristic(NULL);
2367                model->incrementUsed(solver->getColSolution());
2368                cutoff =model->getCutoff();
2369                if (newObjectiveValue >= cutoff)        //  *new* cutoff
2370                  objectiveChange = 1.0e100 ;
2371              }
2372            }
2373          } else if (iStatus==1) {
2374            objectiveChange = 1.0e100 ;
2375          } else {
2376            // Can't say much as we did not finish
2377            choice.finishedDown = false ;
2378          }
2379          choice.downMovement = objectiveChange ;
2380         
2381          // restore bounds
2382          for ( j=0;j<numberColumns;j++) {
2383            if (saveLower[j] != lower[j])
2384              solver->setColLower(j,saveLower[j]);
2385            if (saveUpper[j] != upper[j])
2386              solver->setColUpper(j,saveUpper[j]);
2387          }
2388          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2389          //     choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
2390          //     choice.downMovement,choice.finishedDown,choice.numIntInfeasDown,
2391          //     choice.numObjInfeasDown);
2392         
2393          // repeat the whole exercise, forcing the variable up
2394          decision->saveBranchingObject( choice.possibleBranch);
2395          choice.possibleBranch->branch();
2396          solver->solveFromHotStart() ;
2397          /*
2398            We now have an estimate of objective degradation that we can use for strong
2399            branching. If we're over the cutoff, the variable is monotone up.
2400            If we actually made it to optimality, check for a solution, and if we have
2401            a good one, call setBestSolution to process it. Note that this may reduce the
2402            cutoff, so we check again to see if we can declare this variable monotone.
2403          */
2404          if (solver->isProvenOptimal())
2405            iStatus=0; // optimal
2406          else if (solver->isIterationLimitReached()
2407                   &&!solver->isDualObjectiveLimitReached())
2408            iStatus=2; // unknown
2409          else
2410            iStatus=1; // infeasible
2411          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2412          choice.numItersUp = solver->getIterationCount();
2413          objectiveChange = newObjectiveValue  - objectiveValue_;
2414          decision->updateInformation( solver,this);
2415          if (!iStatus) {
2416            choice.finishedUp = true ;
2417            if (newObjectiveValue>=cutoff) {
2418              objectiveChange = 1.0e100; // say infeasible
2419            } else {
2420              // See if integer solution
2421              if (model->feasibleSolution(choice.numIntInfeasUp,
2422                                          choice.numObjInfeasUp)
2423                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2424                model->setBestSolution(CBC_STRONGSOL,
2425                                       newObjectiveValue,
2426                                       solver->getColSolution()) ;
2427                model->setLastHeuristic(NULL);
2428                model->incrementUsed(solver->getColSolution());
2429                cutoff =model->getCutoff();
2430                if (newObjectiveValue >= cutoff)        //  *new* cutoff
2431                  objectiveChange = 1.0e100 ;
2432              }
2433            }
2434          } else if (iStatus==1) {
2435            objectiveChange = 1.0e100 ;
2436          } else {
2437            // Can't say much as we did not finish
2438            choice.finishedUp = false ;
2439          }
2440          choice.upMovement = objectiveChange ;
2441         
2442          // restore bounds
2443          for ( j=0;j<numberColumns;j++) {
2444            if (saveLower[j] != lower[j])
2445              solver->setColLower(j,saveLower[j]);
2446            if (saveUpper[j] != upper[j])
2447              solver->setColUpper(j,saveUpper[j]);
2448          }
2449         
2450          //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2451          //     choice.objectNumber,iStatus,newObjectiveValue,choice.numItersUp,
2452          //     choice.upMovement,choice.finishedUp,choice.numIntInfeasUp,
2453          //     choice.numObjInfeasUp);
2454        }
2455   
2456        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit2); 
2457        /*
2458          End of evaluation for this candidate variable. Possibilities are:
2459          * Both sides below cutoff; this variable is a candidate for branching.
2460          * Both sides infeasible or above the objective cutoff: no further action
2461          here. Break from the evaluation loop and assume the node will be purged
2462          by the caller.
2463          * One side below cutoff: Install the branch (i.e., fix the variable). Break
2464          from the evaluation loop and assume the node will be reoptimised by the
2465          caller.
2466        */
2467        if (choice.upMovement<1.0e100) {
2468          if(choice.downMovement<1.0e100) {
2469            // In case solution coming in was odd
2470            choice.upMovement = CoinMax(0.0,choice.upMovement);
2471            choice.downMovement = CoinMax(0.0,choice.downMovement);
2472            // feasible - see which best
2473            model->messageHandler()->message(CBC_STRONG,model->messages())
2474              << iObject << iColumn
2475              <<choice.downMovement<<choice.numIntInfeasDown
2476              <<choice.upMovement<<choice.numIntInfeasUp
2477              <<choice.possibleBranch->value()
2478              <<CoinMessageEol;
2479            //if (!stateOfSearch)
2480            //choice.numIntInfeasDown=99999; // temp fudge
2481            int betterWay = decision->betterBranch(choice.possibleBranch,
2482                                                   branch_,
2483                                                   choice.upMovement*changeFactor,
2484                                                   choice.numIntInfeasUp ,
2485                                                   choice.downMovement*changeFactor,
2486                                                   choice.numIntInfeasDown );
2487            if (iDo>=changeStrategy) {
2488              // make less likely
2489              changeStrategy+=numberStrong;
2490              changeFactor *= 0.9;
2491            }
2492            if (betterWay) {
2493              delete branch_;
2494              // C) create branching object
2495              branch_ = choice.possibleBranch;
2496              choice.possibleBranch=NULL;
2497              branch_->way(betterWay);
2498              bestChoice = choice.objectNumber;
2499              whichChoice = iDo;
2500              if (numberStrong<=1) {
2501                delete ws;
2502                ws=NULL;
2503                break;
2504              }
2505            } else {
2506              delete choice.possibleBranch;
2507              if (iDo>=2*numberStrong) {
2508                delete ws;
2509                ws=NULL;
2510                break;
2511              }
2512              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
2513                if (iDo-whichChoice>=numberStrong)
2514                  break; // give up
2515              } else {
2516                if (iDo-whichChoice>=2*numberStrong) {
2517                  delete ws;
2518                  ws=NULL;
2519                  break; // give up
2520                }
2521              }
2522            }
2523          } else {
2524            // up feasible, down infeasible
2525            anyAction=-1;
2526            //printf("Down infeasible for choice %d sequence %d\n",i,
2527            // model->object(choice.objectNumber)->columnNumber());
2528            if (!solveAll) {
2529              choice.possibleBranch->way(1);
2530              choice.possibleBranch->branch();
2531              delete choice.possibleBranch;
2532              delete ws;
2533              ws=NULL;
2534              break;
2535            } else {
2536              choice.fix=1;
2537              fixObject[numberToFix++]=choice;
2538#define FIXNOW
2539#ifdef FIXNOW
2540              double value = ceil(saveSolution[iColumn]);
2541              saveLower[iColumn]=value;
2542              solver->setColLower(iColumn,value);
2543#endif
2544            }
2545          }
2546        } else {
2547          if(choice.downMovement<1.0e100) {
2548            // down feasible, up infeasible
2549            anyAction=-1;
2550            //printf("Up infeasible for choice %d sequence %d\n",i,
2551            // model->object(choice.objectNumber)->columnNumber());
2552            if (!solveAll) {
2553              choice.possibleBranch->way(-1);
2554              choice.possibleBranch->branch();
2555              delete choice.possibleBranch;
2556              delete ws;
2557              ws=NULL;
2558              break;
2559            } else {
2560              choice.fix=-1;
2561              fixObject[numberToFix++]=choice;
2562#ifdef FIXNOW
2563              double value = floor(saveSolution[iColumn]);
2564              saveUpper[iColumn]=value;
2565              solver->setColUpper(iColumn,value);
2566#endif
2567            }
2568          } else {
2569            // neither side feasible
2570            anyAction=-2;
2571            delete choice.possibleBranch;
2572            //printf("Both infeasible for choice %d sequence %d\n",i,
2573            // model->object(choice.objectNumber)->columnNumber());
2574            delete ws;
2575            ws=NULL;
2576            break;
2577          }
2578        }
2579        // Check max time
2580        hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2581                       model->getDblParam(CbcModel::CbcMaximumSeconds));
2582        if (hitMaxTime) {
2583          // make sure rest are fast
2584          doQuickly=true;
2585          for ( int jDo=iDo+1;jDo<numberToDo;jDo++) {
2586            int iObject = whichObject[iDo];
2587            CbcObject * object = model->modifiableObject(iObject);
2588            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2589              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2590            dynamicObject->setNumberBeforeTrust(0);
2591          }
2592          numberTest=0;
2593          distanceToCutoff=-COIN_DBL_MAX;
2594        }
2595      }
2596      if (model->messageHandler()->logLevel()>3||false) { 
2597        if (anyAction==-2)
2598          printf("infeasible\n");
2599        else if(anyAction==-1)
2600          if (!solveAll)
2601            printf("%d fixed\n",numberToFix);
2602          else
2603            printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n",numberToFix,bestChoice,
2604                   iDo,whichChoice,numberToDo);
2605        else
2606          printf("choosing %d iDo %d iChosenWhen %d numberToDo %d\n",bestChoice,
2607                 iDo,whichChoice,numberToDo);
2608      }
2609      if (!skipAll) {
2610        // Delete the snapshot
2611        solver->unmarkHotStart();
2612        // back to normal
2613        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2614        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
2615        // restore basis
2616        solver->setWarmStart(ws);
2617      }
2618      // Unless infeasible we will carry on
2619      // But we could fix anyway
2620      if (numberToFix&&!hitMaxTime) {
2621        if (anyAction==-2) {
2622          // take off
2623          for (i = 0 ; i < numberToFix ; i++) {
2624            delete fixObject[i].possibleBranch;
2625          }
2626        } else {
2627          // apply and take off
2628          for (i = 0 ; i < numberToFix ; i++) {
2629#ifndef FIXNOW
2630            fixObject[i].possibleBranch->way(fixObject[i].fix) ;
2631            fixObject[i].possibleBranch->branch() ;
2632#endif
2633            delete fixObject[i].possibleBranch;
2634          }
2635          bool feasible=true;
2636#if ACTION <2
2637          // can do quick optimality check
2638          int easy=2;
2639          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
2640          solver->resolve() ;
2641          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2642          feasible = solver->isProvenOptimal();
2643          if (feasible) {
2644            anyAction=0;
2645            memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2646            model->reserveCurrentSolution(saveSolution);
2647            memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
2648            memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
2649            // See if candidate still possible
2650            if (branch_) {
2651              const CbcObject * object = model->object(bestChoice);
2652              int preferredWay;
2653              double infeasibility = object->infeasibility(preferredWay);
2654              if (!infeasibility) {
2655                // take out
2656                delete branch_;
2657                branch_=NULL;
2658              } else {
2659                branch_->way(preferredWay);
2660              }
2661            }
2662          } else {
2663            anyAction=-2;
2664            finished=true;
2665          }
2666#endif
2667          // If  fixed then round again
2668          if (!branch_&&anyAction!=-2) {
2669            finished=false;
2670          }
2671          // If these in then different action
2672#if ACTION == 1
2673          if (!anyAction)
2674            anyAction=-1;
2675          finished=true;
2676#endif
2677        }
2678      }
2679      delete ws;
2680    }
2681  }
2682 
2683  // Set guessed solution value
2684  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
2685/*
2686  Cleanup, then we're finished
2687*/
2688  if (!model->branchingMethod())
2689    delete decision;
2690   
2691  delete [] fixObject;
2692  delete [] sort;
2693  delete [] whichObject;
2694  delete [] objectMark;
2695  delete [] saveLower;
2696  delete [] saveUpper;
2697  delete [] upEstimate;
2698  delete [] downEstimate;
2699 
2700  // restore solution
2701  solver->setColSolution(saveSolution);
2702  model->reserveCurrentSolution(saveSolution);
2703  delete [] saveSolution;
2704  return anyAction;
2705}
2706
2707
2708CbcNode::CbcNode(const CbcNode & rhs) 
2709{ 
2710#ifdef CHECK_NODE
2711  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
2712#endif
2713  if (rhs.nodeInfo_)
2714    nodeInfo_ = rhs.nodeInfo_->clone();
2715  else
2716    nodeInfo_=NULL;
2717  objectiveValue_=rhs.objectiveValue_;
2718  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
2719  if (rhs.branch_)
2720    branch_=rhs.branch_->clone();
2721  else
2722    branch_=NULL;
2723  depth_ = rhs.depth_;
2724  numberUnsatisfied_ = rhs.numberUnsatisfied_;
2725}
2726
2727CbcNode &
2728CbcNode::operator=(const CbcNode & rhs)
2729{
2730  if (this != &rhs) {
2731    delete nodeInfo_;
2732    if (nodeInfo_)
2733      nodeInfo_ = rhs.nodeInfo_->clone();
2734    else
2735      nodeInfo_ = NULL;
2736    objectiveValue_=rhs.objectiveValue_;
2737    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
2738    if (rhs.branch_)
2739      branch_=rhs.branch_->clone();
2740    else
2741      branch_=NULL,
2742    depth_ = rhs.depth_;
2743    numberUnsatisfied_ = rhs.numberUnsatisfied_;
2744  }
2745  return *this;
2746}
2747
2748
2749CbcNode::~CbcNode ()
2750{
2751#ifdef CHECK_NODE
2752  if (nodeInfo_) 
2753    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
2754         this,nodeInfo_,nodeInfo_->numberPointingToThis());
2755  else
2756    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
2757         this,nodeInfo_);
2758#endif
2759  if (nodeInfo_) {
2760    nodeInfo_->nullOwner();
2761    int numberToDelete=nodeInfo_->numberBranchesLeft();
2762    //    CbcNodeInfo * parent = nodeInfo_->parent();
2763    if (nodeInfo_->decrement(numberToDelete)==0) {
2764      delete nodeInfo_;
2765    } else {
2766      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent);
2767      // anyway decrement parent
2768      //if (parent)
2769      ///parent->decrement(1);
2770    }
2771  }
2772  delete branch_;
2773}
2774// Decrement  active cut counts
2775void 
2776CbcNode::decrementCuts(int change)
2777{
2778  if(nodeInfo_) {
2779    nodeInfo_->decrementCuts(change);
2780  }
2781}
2782void 
2783CbcNode::decrementParentCuts(int change)
2784{
2785  if(nodeInfo_) {
2786    nodeInfo_->decrementParentCuts(change);
2787  }
2788}
2789
2790/*
2791  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
2792  in the attached nodeInfo_.
2793*/
2794void
2795CbcNode::initializeInfo()
2796{
2797  assert(nodeInfo_ && branch_) ;
2798  nodeInfo_->initializeInfo(branch_->numberBranches());
2799}
2800// Nulls out node info
2801void 
2802CbcNode::nullNodeInfo()
2803{
2804  nodeInfo_=NULL;
2805}
2806
2807int
2808CbcNode::branch()
2809{
2810  double changeInGuessed=branch_->branch(true);
2811  guessedObjectiveValue_+= changeInGuessed;
2812  //#define PRINTIT
2813#ifdef PRINTIT
2814  int numberLeft = nodeInfo_->numberBranchesLeft();
2815  CbcNodeInfo * parent = nodeInfo_->parent();
2816  int parentNodeNumber = -1;
2817  //CbcBranchingObject * object1 = branch_->object_;
2818  //CbcObject * object = object1->
2819  //int sequence = object->columnNumber);
2820  int id=-1;
2821  double value=0.0;
2822  if (branch_) {
2823    id = branch_->variable();
2824    value = branch_->value();
2825  }
2826  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
2827  if (parent)
2828    parentNodeNumber = parent->nodeNumber();
2829  printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
2830         nodeInfo_->nodeNumber(),(numberLeft==2) ? "leftBranch" : "rightBranch",
2831         way(),depth_,parentNodeNumber);
2832#endif
2833  return nodeInfo_->branchedOn();
2834}
Note: See TracBrowser for help on using the repository browser.