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

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

for nonlinear

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