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

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

try and make faster

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