source: stable/2.2/Cbc/src/CbcBranchActual.cpp @ 1050

Last change on this file since 1050 was 1050, checked in by forrest, 11 years ago

sos bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 142.5 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cstdlib>
9#include <cmath>
10#include <cfloat>
11//#define CBC_DEBUG
12
13#include "CoinTypes.hpp"
14#include "OsiSolverInterface.hpp"
15#include "OsiSolverBranch.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcBranchActual.hpp"
19#include "CoinSort.hpp"
20#include "CoinError.hpp"
21
22//##############################################################################
23
24// Default Constructor
25CbcClique::CbcClique ()
26  : CbcObject(),
27    numberMembers_(0),
28    numberNonSOSMembers_(0),
29    members_(NULL),
30    type_(NULL),
31    cliqueType_(-1),
32    slack_(-1)
33{
34}
35
36// Useful constructor (which are integer indices)
37CbcClique::CbcClique (CbcModel * model, int cliqueType, int numberMembers,
38           const int * which, const char * type, int identifier,int slack)
39  : CbcObject(model)
40{
41  id_=identifier;
42  numberMembers_=numberMembers;
43  if (numberMembers_) {
44    members_ = new int[numberMembers_];
45    memcpy(members_,which,numberMembers_*sizeof(int));
46    type_ = new char[numberMembers_];
47    if (type) {
48      memcpy(type_,type,numberMembers_*sizeof(char));
49    } else {
50      for (int i=0;i<numberMembers_;i++)
51        type_[i]=1;
52    }
53  } else {
54    members_ = NULL;
55    type_ = NULL;
56  }
57  // Find out how many non sos
58  int i;
59  numberNonSOSMembers_=0;
60  for (i=0;i<numberMembers_;i++)
61    if (!type_[i])
62      numberNonSOSMembers_++;
63  cliqueType_ = cliqueType;
64  slack_ = slack;
65}
66
67// Copy constructor
68CbcClique::CbcClique ( const CbcClique & rhs)
69  :CbcObject(rhs)
70{
71  numberMembers_ = rhs.numberMembers_;
72  numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
73  if (numberMembers_) {
74    members_ = new int[numberMembers_];
75    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
76    type_ = new char[numberMembers_];
77    memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
78  } else {
79    members_ = NULL;
80    type_ = NULL;
81  }
82  cliqueType_ = rhs.cliqueType_;
83  slack_ = rhs.slack_;
84}
85
86// Clone
87CbcObject *
88CbcClique::clone() const
89{
90  return new CbcClique(*this);
91}
92
93// Assignment operator
94CbcClique & 
95CbcClique::operator=( const CbcClique& rhs)
96{
97  if (this!=&rhs) {
98    CbcObject::operator=(rhs);
99    delete [] members_;
100    delete [] type_;
101    numberMembers_ = rhs.numberMembers_;
102    numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
103    if (numberMembers_) {
104      members_ = new int[numberMembers_];
105      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
106      type_ = new char[numberMembers_];
107      memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
108    } else {
109      members_ = NULL;
110      type_ = NULL;
111    }
112    cliqueType_ = rhs.cliqueType_;
113    slack_ = rhs.slack_;
114  }
115  return *this;
116}
117
118// Destructor
119CbcClique::~CbcClique ()
120{
121  delete [] members_;
122  delete [] type_;
123}
124
125// Infeasibility - large is 0.5
126double 
127CbcClique::infeasibility(int & preferredWay) const
128{
129  int numberUnsatis=0, numberFree=0;
130  int j;
131  const int * integer = model_->integerVariable();
132  OsiSolverInterface * solver = model_->solver();
133  const double * solution = model_->testSolution();
134  const double * lower = solver->getColLower();
135  const double * upper = solver->getColUpper();
136  double largestValue=0.0;
137  double integerTolerance = 
138    model_->getDblParam(CbcModel::CbcIntegerTolerance);
139  double * sort = new double[numberMembers_];
140
141  double slackValue=0.0;
142  for (j=0;j<numberMembers_;j++) {
143    int sequence = members_[j];
144    int iColumn = integer[sequence];
145    double value = solution[iColumn];
146    value = CoinMax(value, lower[iColumn]);
147    value = CoinMin(value, upper[iColumn]);
148    double nearest = floor(value+0.5);
149    double distance = fabs(value-nearest);
150    if (distance>integerTolerance) {
151      if (!type_[j])
152        value = 1.0-value; // non SOS
153      // if slack then choose that
154      if (j==slack_&&value>0.05)
155        slackValue = value;
156      largestValue = CoinMax(value,largestValue);
157      sort[numberUnsatis++]=-value;
158    } else if (upper[iColumn]>lower[iColumn]) {
159      numberFree++;
160    }
161  }
162  preferredWay=1;
163  double otherWay = 0.0;
164  if (numberUnsatis) {
165    // sort
166    std::sort(sort,sort+numberUnsatis);
167    for (j=0;j<numberUnsatis;j++) {
168      if ((j&1)!=0)
169        otherWay += -sort[j];
170    }
171    // Need to think more
172    double value = 0.2*numberUnsatis+0.01*(numberMembers_-numberFree);
173    if (fabs(largestValue-0.5)<0.1) {
174      // close to half
175      value +=0.1;
176    }
177    if (slackValue) {
178      // branching on slack
179      value += slackValue;
180    }
181    // scale other way
182    otherWay *= value/(1.0-otherWay);
183    delete [] sort;
184    return value;
185  } else {
186    delete [] sort;
187    return 0.0; // satisfied
188  }
189}
190
191// This looks at solution and sets bounds to contain solution
192void 
193CbcClique::feasibleRegion()
194{
195  int j;
196  const int * integer = model_->integerVariable();
197  OsiSolverInterface * solver = model_->solver();
198  const double * solution = model_->testSolution();
199  const double * lower = solver->getColLower();
200  const double * upper = solver->getColUpper();
201#ifndef NDEBUG
202  double integerTolerance = 
203    model_->getDblParam(CbcModel::CbcIntegerTolerance);
204#endif 
205  for (j=0;j<numberMembers_;j++) {
206    int sequence = members_[j];
207    int iColumn = integer[sequence];
208    double value = solution[iColumn];
209    value = CoinMax(value, lower[iColumn]);
210    value = CoinMin(value, upper[iColumn]);
211    double nearest = floor(value+0.5);
212#ifndef NDEBUG
213    double distance = fabs(value-nearest);
214    assert(distance<=integerTolerance);
215#endif
216    solver->setColLower(iColumn,nearest);
217    solver->setColUpper(iColumn,nearest);
218  }
219}
220// Redoes data when sequence numbers change
221void 
222CbcClique::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
223{
224  model_=model;
225  int n2=0;
226  for (int j=0;j<numberMembers_;j++) {
227    int iColumn = members_[j];
228    int i;
229    for (i=0;i<numberColumns;i++) {
230      if (originalColumns[i]==iColumn)
231        break;
232    }
233    if (i<numberColumns) {
234      members_[n2]=i;
235      type_[n2++]=type_[j];
236    }
237  }
238  if (n2<numberMembers_) {
239    //printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
240    numberMembers_=n2;
241  }
242  // Find out how many non sos
243  int i;
244  numberNonSOSMembers_=0;
245  for (i=0;i<numberMembers_;i++)
246    if (!type_[i])
247      numberNonSOSMembers_++;
248}
249
250
251// Creates a branching object
252CbcBranchingObject * 
253CbcClique::createBranch(int way) 
254{
255  int numberUnsatis=0;
256  int j;
257  int nUp=0;
258  int nDown=0;
259  int numberFree=numberMembers_;
260  const int * integer = model_->integerVariable();
261  OsiSolverInterface * solver = model_->solver();
262  const double * solution = model_->testSolution();
263  const double * lower = solver->getColLower();
264  const double * upper = solver->getColUpper();
265  int * upList = new int[numberMembers_];
266  int * downList = new int[numberMembers_];
267  double * sort = new double[numberMembers_];
268  double integerTolerance = 
269      model_->getDblParam(CbcModel::CbcIntegerTolerance);
270
271  double slackValue=0.0;
272  for (j=0;j<numberMembers_;j++) {
273    int sequence = members_[j];
274    int iColumn = integer[sequence];
275    double value = solution[iColumn];
276    value = CoinMax(value, lower[iColumn]);
277    value = CoinMin(value, upper[iColumn]);
278    double nearest = floor(value+0.5);
279    double distance = fabs(value-nearest);
280    if (distance>integerTolerance) {
281      if (!type_[j])
282        value = 1.0-value; // non SOS
283      // if slack then choose that
284      if (j==slack_&&value>0.05)
285        slackValue = value;
286      value = -value; // for sort
287      upList[numberUnsatis]=j;
288      sort[numberUnsatis++]=value;
289    } else if (upper[iColumn]>lower[iColumn]) {
290      upList[--numberFree]=j;
291    }
292  }
293  assert (numberUnsatis);
294  if (!slackValue) {
295    // sort
296    CoinSort_2(sort,sort+numberUnsatis,upList);
297    // put first in up etc
298    int kWay=1;
299    for (j=0;j<numberUnsatis;j++) {
300      if (kWay>0) 
301        upList[nUp++]=upList[j];
302      else
303        downList[nDown++]=upList[j];
304      kWay = -kWay;
305    }
306    for (j=numberFree;j<numberMembers_;j++) {
307      if (kWay>0)
308        upList[nUp++]=upList[j];
309      else
310        downList[nDown++]=upList[j];
311      kWay = -kWay;
312    }
313  } else {
314    // put slack to 0 in first way
315    nUp = 1;
316    upList[0]=slack_;
317    for (j=0;j<numberUnsatis;j++) {
318      downList[nDown++]=upList[j];
319    }
320    for (j=numberFree;j<numberMembers_;j++) {
321      downList[nDown++]=upList[j];
322    }
323  }
324  // create object
325  CbcBranchingObject * branch;
326  if (numberMembers_ <=64)
327     branch = new CbcCliqueBranchingObject(model_,this,way,
328                                         nDown,downList,nUp,upList);
329  else
330    branch = new CbcLongCliqueBranchingObject(model_,this,way,
331                                            nDown,downList,nUp,upList);
332  delete [] upList;
333  delete [] downList;
334  delete [] sort;
335  return branch;
336}
337
338//##############################################################################
339
340// Default Constructor
341CbcSOS::CbcSOS ()
342  : CbcObject(),
343    members_(NULL),
344    weights_(NULL),
345    shadowEstimateDown_(1.0),
346    shadowEstimateUp_(1.0),
347    downDynamicPseudoRatio_(0.0),
348    upDynamicPseudoRatio_(0.0),
349    numberTimesDown_(0),
350    numberTimesUp_(0),
351    numberMembers_(0),
352    sosType_(-1),
353    integerValued_(false)
354{
355}
356
357// Useful constructor (which are indices)
358CbcSOS::CbcSOS (CbcModel * model,  int numberMembers,
359           const int * which, const double * weights, int identifier,int type)
360  : CbcObject(model),
361    shadowEstimateDown_(1.0),
362    shadowEstimateUp_(1.0),
363    downDynamicPseudoRatio_(0.0),
364    upDynamicPseudoRatio_(0.0),
365    numberTimesDown_(0),
366    numberTimesUp_(0),
367    numberMembers_(numberMembers),
368    sosType_(type)
369{
370  id_=identifier;
371  integerValued_ = type==1;
372  if (integerValued_) {
373    // check all members integer
374    OsiSolverInterface * solver = model->solver();
375    if (solver) {
376      for (int i=0;i<numberMembers_;i++) {
377        if (!solver->isInteger(which[i]))
378          integerValued_=false;
379      }
380    } else {
381      // can't tell
382      integerValued_=false;
383    }
384  }
385  if (numberMembers_) {
386    members_ = new int[numberMembers_];
387    weights_ = new double[numberMembers_];
388    memcpy(members_,which,numberMembers_*sizeof(int));
389    if (weights) {
390      memcpy(weights_,weights,numberMembers_*sizeof(double));
391    } else {
392      for (int i=0;i<numberMembers_;i++)
393        weights_[i]=i;
394    }
395    // sort so weights increasing
396    CoinSort_2(weights_,weights_+numberMembers_,members_);
397    double last = -COIN_DBL_MAX;
398    int i;
399    for (i=0;i<numberMembers_;i++) {
400      double possible = CoinMax(last+1.0e-10,weights_[i]);
401      weights_[i] = possible;
402      last=possible;
403    }
404  } else {
405    members_ = NULL;
406    weights_ = NULL;
407  }
408  assert (sosType_>0&&sosType_<3);
409}
410
411// Copy constructor
412CbcSOS::CbcSOS ( const CbcSOS & rhs)
413  :CbcObject(rhs)
414{
415  shadowEstimateDown_ = rhs.shadowEstimateDown_;
416  shadowEstimateUp_ = rhs.shadowEstimateUp_;
417  downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
418  upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
419  numberTimesDown_ = rhs.numberTimesDown_;
420  numberTimesUp_ = rhs.numberTimesUp_;
421  numberMembers_ = rhs.numberMembers_;
422  sosType_ = rhs.sosType_;
423  integerValued_ = rhs.integerValued_;
424  if (numberMembers_) {
425    members_ = new int[numberMembers_];
426    weights_ = new double[numberMembers_];
427    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
428    memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
429  } else {
430    members_ = NULL;
431    weights_ = NULL;
432  }
433}
434
435// Clone
436CbcObject *
437CbcSOS::clone() const
438{
439  return new CbcSOS(*this);
440}
441
442// Assignment operator
443CbcSOS & 
444CbcSOS::operator=( const CbcSOS& rhs)
445{
446  if (this!=&rhs) {
447    CbcObject::operator=(rhs);
448    delete [] members_;
449    delete [] weights_;
450    shadowEstimateDown_ = rhs.shadowEstimateDown_;
451    shadowEstimateUp_ = rhs.shadowEstimateUp_;
452    downDynamicPseudoRatio_ = rhs.downDynamicPseudoRatio_;
453    upDynamicPseudoRatio_ = rhs.upDynamicPseudoRatio_;
454    numberTimesDown_ = rhs.numberTimesDown_;
455    numberTimesUp_ = rhs.numberTimesUp_;
456    numberMembers_ = rhs.numberMembers_;
457    sosType_ = rhs.sosType_;
458    integerValued_ = rhs.integerValued_;
459    if (numberMembers_) {
460      members_ = new int[numberMembers_];
461      weights_ = new double[numberMembers_];
462      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
463      memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
464    } else {
465      members_ = NULL;
466      weights_ = NULL;
467    }
468  }
469  return *this;
470}
471
472// Destructor
473CbcSOS::~CbcSOS ()
474{
475  delete [] members_;
476  delete [] weights_;
477}
478
479// Infeasibility - large is 0.5
480double 
481CbcSOS::infeasibility(int & preferredWay) const
482{
483  int j;
484  int firstNonZero=-1;
485  int lastNonZero = -1;
486  OsiSolverInterface * solver = model_->solver();
487  const double * solution = model_->testSolution();
488  //const double * lower = solver->getColLower();
489  const double * upper = solver->getColUpper();
490  //double largestValue=0.0;
491  double integerTolerance = 
492    model_->getDblParam(CbcModel::CbcIntegerTolerance);
493  double weight = 0.0;
494  double sum =0.0;
495
496  // check bounds etc
497  double lastWeight=-1.0e100;
498  for (j=0;j<numberMembers_;j++) {
499    int iColumn = members_[j];
500    if (lastWeight>=weights_[j]-1.0e-7)
501      throw CoinError("Weights too close together in SOS","infeasibility","CbcSOS");
502    double value = CoinMax(0.0,solution[iColumn]);
503    sum += value;
504    if (value>integerTolerance&&upper[iColumn]) {
505      // Possibly due to scaling a fixed variable might slip through
506      if (value>upper[iColumn]) {
507        value=upper[iColumn];
508        // Could change to #ifdef CBC_DEBUG
509#ifndef NDEBUG
510        if (model_->messageHandler()->logLevel()>2)
511          printf("** Variable %d (%d) has value %g and upper bound of %g\n",
512                 iColumn,j,value,upper[iColumn]);
513#endif
514      } 
515      weight += weights_[j]*value;
516      if (firstNonZero<0)
517        firstNonZero=j;
518      lastNonZero=j;
519    }
520  }
521  preferredWay=1;
522  if (lastNonZero-firstNonZero>=sosType_) {
523    // find where to branch
524    assert (sum>0.0);
525    weight /= sum;
526    //int iWhere;
527    //for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++)
528    //if (weight<weights_[iWhere+1])
529    //break;
530    // probably best to use pseudo duals
531    double value = lastNonZero-firstNonZero+1;
532    value *= 0.5/((double) numberMembers_);
533    // adjust??
534    return value;
535  } else {
536    return 0.0; // satisfied
537  }
538}
539// Infeasibility - large is 0.5
540double 
541CbcSOS::infeasibility(const OsiBranchingInformation * info, 
542                      int & preferredWay) const
543{
544  int j;
545  int firstNonZero=-1;
546  int lastNonZero = -1;
547  OsiSolverInterface * solver = model_->solver();
548  const double * solution = model_->testSolution();
549  //const double * lower = solver->getColLower();
550  const double * upper = solver->getColUpper();
551  //double largestValue=0.0;
552  double integerTolerance = 
553    model_->getDblParam(CbcModel::CbcIntegerTolerance);
554  double weight = 0.0;
555  double sum =0.0;
556
557  // check bounds etc
558  double lastWeight=-1.0e100;
559  for (j=0;j<numberMembers_;j++) {
560    int iColumn = members_[j];
561    if (lastWeight>=weights_[j]-1.0e-7)
562      throw CoinError("Weights too close together in SOS","infeasibility","CbcSOS");
563    double value = CoinMax(0.0,solution[iColumn]);
564    sum += value;
565    if (value>integerTolerance&&upper[iColumn]) {
566      // Possibly due to scaling a fixed variable might slip through
567      if (value>upper[iColumn]) {
568        value=upper[iColumn];
569        // Could change to #ifdef CBC_DEBUG
570#ifndef NDEBUG
571        if (model_->messageHandler()->logLevel()>2)
572          printf("** Variable %d (%d) has value %g and upper bound of %g\n",
573                 iColumn,j,value,upper[iColumn]);
574#endif
575      } 
576      weight += weights_[j]*value;
577      if (firstNonZero<0)
578        firstNonZero=j;
579      lastNonZero=j;
580    }
581  }
582  preferredWay=1;
583  if (lastNonZero-firstNonZero>=sosType_) {
584    // find where to branch
585    assert (sum>0.0);
586    weight /= sum;
587    if (info->defaultDual_>=0.0&&info->usefulRegion_) {
588      assert (sosType_==1);
589      int iWhere;
590      for (iWhere=firstNonZero;iWhere<lastNonZero-1;iWhere++) {
591        if (weight<weights_[iWhere+1]) {
592          break;
593        }
594      }
595      int jColumnDown = members_[iWhere];
596      int jColumnUp = members_[iWhere+1];
597      int n=0;
598      CoinBigIndex j;
599      double objMove = info->objective_[jColumnDown];
600      for (j=info->columnStart_[jColumnDown];
601           j<info->columnStart_[jColumnDown]+info->columnLength_[jColumnDown];j++) {
602        double value = info->elementByColumn_[j];
603        int iRow = info->row_[j];
604        info->indexRegion_[n++]=iRow;
605        info->usefulRegion_[iRow]=value;
606      }
607      for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) {
608        int jColumn = members_[iWhere];
609        double solValue = info->solution_[jColumn];
610        if (!solValue)
611          continue;
612        objMove -= info->objective_[jColumn]*solValue;
613        for (j=info->columnStart_[jColumn];
614             j<info->columnStart_[jColumn]+info->columnLength_[jColumn];j++) {
615          double value = -info->elementByColumn_[j]*solValue;
616          int iRow = info->row_[j];
617          double oldValue = info->usefulRegion_[iRow];
618          if (!oldValue) {
619            info->indexRegion_[n++]=iRow;
620          } else {
621            value += oldValue;
622            if (!value)
623              value = 1.0e-100;
624          }
625          info->usefulRegion_[iRow]=value;
626        }
627      }
628      const double * pi = info->pi_;
629      const double * activity = info->rowActivity_;
630      const double * lower = info->rowLower_;
631      const double * upper = info->rowUpper_;
632      double tolerance = info->primalTolerance_;
633      double direction = info->direction_;
634      shadowEstimateDown_ = objMove*direction;
635      bool infeasible=false;
636      for (int k=0;k<n;k++) {
637        int iRow = info->indexRegion_[k];
638        double movement=info->usefulRegion_[iRow];
639        // not this time info->usefulRegion_[iRow]=0.0;
640        if (lower[iRow]<-1.0e20) 
641          assert (pi[iRow]<=1.0e-3);
642        if (upper[iRow]>1.0e20) 
643          assert (pi[iRow]>=-1.0e-3);
644        double valueP = pi[iRow]*direction;
645        // if move makes infeasible then make at least default
646        double newValue = activity[iRow] + movement;
647        if (newValue>upper[iRow]+tolerance||newValue<lower[iRow]-tolerance) {
648          shadowEstimateDown_ += fabs(movement)*CoinMax(fabs(valueP),info->defaultDual_);
649          infeasible=true;
650        }
651      }
652      if (shadowEstimateDown_<info->integerTolerance_) {
653        if (!infeasible) {
654          shadowEstimateDown_=1.0e-10;
655#ifdef COIN_DEVELOP
656          printf("zero pseudoShadowPrice\n");
657#endif
658        } else
659          shadowEstimateDown_ = info->integerTolerance_;
660      }
661      // And other way
662      // take off
663      objMove -= info->objective_[jColumnDown];
664      for (j=info->columnStart_[jColumnDown];
665           j<info->columnStart_[jColumnDown]+info->columnLength_[jColumnDown];j++) {
666        double value = -info->elementByColumn_[j];
667        int iRow = info->row_[j];
668        double oldValue = info->usefulRegion_[iRow];
669        if (!oldValue) {
670          info->indexRegion_[n++]=iRow;
671          } else {
672          value += oldValue;
673          if (!value)
674            value = 1.0e-100;
675        }
676        info->usefulRegion_[iRow]=value;
677      }
678      // add on
679      objMove += info->objective_[jColumnUp];
680      for (j=info->columnStart_[jColumnUp];
681           j<info->columnStart_[jColumnUp]+info->columnLength_[jColumnUp];j++) {
682        double value = info->elementByColumn_[j];
683        int iRow = info->row_[j];
684        double oldValue = info->usefulRegion_[iRow];
685        if (!oldValue) {
686          info->indexRegion_[n++]=iRow;
687          } else {
688          value += oldValue;
689          if (!value)
690            value = 1.0e-100;
691        }
692        info->usefulRegion_[iRow]=value;
693      }
694      shadowEstimateUp_ = objMove*direction;
695      infeasible=false;
696      for (int k=0;k<n;k++) {
697        int iRow = info->indexRegion_[k];
698        double movement=info->usefulRegion_[iRow];
699        info->usefulRegion_[iRow]=0.0;
700        if (lower[iRow]<-1.0e20) 
701          assert (pi[iRow]<=1.0e-3);
702        if (upper[iRow]>1.0e20) 
703          assert (pi[iRow]>=-1.0e-3);
704        double valueP = pi[iRow]*direction;
705        // if move makes infeasible then make at least default
706        double newValue = activity[iRow] + movement;
707        if (newValue>upper[iRow]+tolerance||newValue<lower[iRow]-tolerance) {
708          shadowEstimateUp_ += fabs(movement)*CoinMax(fabs(valueP),info->defaultDual_);
709          infeasible=true;
710        }
711      }
712      if (shadowEstimateUp_<info->integerTolerance_) {
713        if (!infeasible) {
714          shadowEstimateUp_=1.0e-10;
715#ifdef COIN_DEVELOP
716          printf("zero pseudoShadowPrice\n");
717#endif
718        } else
719          shadowEstimateUp_ = info->integerTolerance_;
720      }
721      // adjust
722      double downCost = shadowEstimateDown_;
723      double upCost = shadowEstimateUp_;
724      if (numberTimesDown_) 
725        downCost *= downDynamicPseudoRatio_/
726          ((double) numberTimesDown_);
727      if (numberTimesUp_) 
728        upCost *= upDynamicPseudoRatio_/
729          ((double) numberTimesUp_);
730#define WEIGHT_AFTER 0.7
731#define WEIGHT_BEFORE 0.1
732      int stateOfSearch = model_->stateOfSearch()%10;
733      double returnValue=0.0;
734      double minValue = CoinMin(downCost,upCost);
735      double maxValue = CoinMax(downCost,upCost);
736      if (stateOfSearch<=2) {
737        // no branching solution
738        returnValue = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
739      } else {
740        returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
741      }
742#ifdef PRINT_SHADOW
743      printf("%d id - down %d %g up %d %g shadow %g, %g returned %g\n",
744             id_,numberTimesDown_,downDynamicPseudoRatio_,
745             numberTimesUp_,upDynamicPseudoRatio_,shadowEstimateDown_,
746             shadowEstimateUp_,returnValue);
747#endif
748      return returnValue;
749    } else {
750      double value = lastNonZero-firstNonZero+1;
751      value *= 0.5/((double) numberMembers_);
752      return value;
753    }
754  } else {
755    return 0.0; // satisfied
756  }
757}
758
759// This looks at solution and sets bounds to contain solution
760void 
761CbcSOS::feasibleRegion()
762{
763  int j;
764  int firstNonZero=-1;
765  int lastNonZero = -1;
766  OsiSolverInterface * solver = model_->solver();
767  const double * solution = model_->testSolution();
768  //const double * lower = solver->getColLower();
769  const double * upper = solver->getColUpper();
770  double integerTolerance = 
771    model_->getDblParam(CbcModel::CbcIntegerTolerance);
772  double weight = 0.0;
773  double sum =0.0;
774
775  for (j=0;j<numberMembers_;j++) {
776    int iColumn = members_[j];
777    double value = CoinMax(0.0,solution[iColumn]);
778    sum += value;
779    if (value>integerTolerance&&upper[iColumn]) {
780      weight += weights_[j]*value;
781      if (firstNonZero<0)
782        firstNonZero=j;
783      lastNonZero=j;
784    }
785  }
786  assert (lastNonZero-firstNonZero<sosType_) ;
787  for (j=0;j<firstNonZero;j++) {
788    int iColumn = members_[j];
789    solver->setColUpper(iColumn,0.0);
790  }
791  for (j=lastNonZero+1;j<numberMembers_;j++) {
792    int iColumn = members_[j];
793    solver->setColUpper(iColumn,0.0);
794  }
795}
796// Redoes data when sequence numbers change
797void 
798CbcSOS::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
799{
800  model_=model;
801  int n2=0;
802  for (int j=0;j<numberMembers_;j++) {
803    int iColumn = members_[j];
804    int i;
805    for (i=0;i<numberColumns;i++) {
806      if (originalColumns[i]==iColumn)
807        break;
808    }
809    if (i<numberColumns) {
810      members_[n2]=i;
811      weights_[n2++]=weights_[j];
812    }
813  }
814  if (n2<numberMembers_) {
815    //printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
816    numberMembers_=n2;
817  }
818}
819
820
821// Creates a branching object
822CbcBranchingObject * 
823CbcSOS::createBranch(int way) 
824{
825  int j;
826  const double * solution = model_->testSolution();
827  double integerTolerance = 
828      model_->getDblParam(CbcModel::CbcIntegerTolerance);
829  OsiSolverInterface * solver = model_->solver();
830  const double * upper = solver->getColUpper();
831  int firstNonFixed=-1;
832  int lastNonFixed=-1;
833  int firstNonZero=-1;
834  int lastNonZero = -1;
835  double weight = 0.0;
836  double sum =0.0;
837  for (j=0;j<numberMembers_;j++) {
838    int iColumn = members_[j];
839    if (upper[iColumn]) {
840      double value = CoinMax(0.0,solution[iColumn]);
841      sum += value;
842      if (firstNonFixed<0)
843        firstNonFixed=j;
844      lastNonFixed=j;
845      if (value>integerTolerance) {
846        weight += weights_[j]*value;
847        if (firstNonZero<0)
848          firstNonZero=j;
849        lastNonZero=j;
850      }
851    }
852  }
853  assert (lastNonZero-firstNonZero>=sosType_) ;
854  // find where to branch
855  assert (sum>0.0);
856  weight /= sum;
857  int iWhere;
858  double separator=0.0;
859  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
860    if (weight<weights_[iWhere+1])
861      break;
862  if (sosType_==1) {
863    // SOS 1
864    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
865  } else {
866    // SOS 2
867    if (iWhere==firstNonFixed)
868      iWhere++;;
869    if (iWhere==lastNonFixed-1)
870      iWhere = lastNonFixed-2;
871    separator = weights_[iWhere+1];
872  }
873  // create object
874  CbcBranchingObject * branch;
875  branch = new CbcSOSBranchingObject(model_,this,way,separator);
876  branch->setOriginalObject(this);
877  return branch;
878}
879/* Pass in information on branch just done and create CbcObjectUpdateData instance.
880   If object does not need data then backward pointer will be NULL.
881   Assumes can get information from solver */
882CbcObjectUpdateData
883CbcSOS::createUpdateInformation(const OsiSolverInterface * solver, 
884                                const CbcNode * node,
885                                const CbcBranchingObject * branchingObject)
886{
887  double originalValue=node->objectiveValue();
888  int originalUnsatisfied = node->numberUnsatisfied();
889  double objectiveValue = solver->getObjValue()*solver->getObjSense();
890  int unsatisfied=0;
891  int i;
892  //might be base model - doesn't matter
893  int numberIntegers = model_->numberIntegers();;
894  const double * solution = solver->getColSolution();
895  //const double * lower = solver->getColLower();
896  //const double * upper = solver->getColUpper();
897  double change = CoinMax(0.0,objectiveValue-originalValue);
898  int iStatus;
899  if (solver->isProvenOptimal())
900    iStatus=0; // optimal
901  else if (solver->isIterationLimitReached()
902           &&!solver->isDualObjectiveLimitReached())
903    iStatus=2; // unknown
904  else
905    iStatus=1; // infeasible
906
907  bool feasible = iStatus!=1;
908  if (feasible) {
909    double integerTolerance = 
910      model_->getDblParam(CbcModel::CbcIntegerTolerance);
911    const int * integerVariable = model_->integerVariable();
912    for (i=0;i<numberIntegers;i++) {
913      int j=integerVariable[i];
914      double value = solution[j];
915      double nearest = floor(value+0.5);
916      if (fabs(value-nearest)>integerTolerance) 
917        unsatisfied++;
918    }
919  }
920  int way = branchingObject->way();
921  way = - way; // because after branch so moved on
922  double value = branchingObject->value();
923  CbcObjectUpdateData newData (this, way,
924                               change, iStatus,
925                               originalUnsatisfied-unsatisfied,value);
926  newData.originalObjective_ = originalValue;
927  // Solvers know about direction
928  double direction = solver->getObjSense();
929  solver->getDblParam(OsiDualObjectiveLimit,newData.cutoff_);
930  newData.cutoff_ *= direction;
931  return newData;
932}
933// Update object by CbcObjectUpdateData
934void 
935CbcSOS::updateInformation(const CbcObjectUpdateData & data)
936{
937  bool feasible = data.status_!=1;
938  int way = data.way_;
939  //double value = data.branchingValue_;
940  double originalValue = data.originalObjective_;
941  double change = data.change_;
942  if (way<0) {
943    // down
944    if (!feasible) {
945      double distanceToCutoff=0.0;
946      //double objectiveValue = model_->getCurrentMinimizationObjValue();
947      distanceToCutoff =  model_->getCutoff()  - originalValue;
948      if (distanceToCutoff<1.0e20) 
949        change = distanceToCutoff*2.0;
950      else 
951        change = (downDynamicPseudoRatio_*shadowEstimateDown_+1.0e-3)*10.0;
952    } 
953    change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
954#ifdef PRINT_SHADOW
955    if (numberTimesDown_)
956      printf("Updating id %d - down change %g (true %g) - ndown %d estimated change %g - raw shadow estimate %g\n",
957             id_,change,data.change_,numberTimesDown_,shadowEstimateDown_*
958             (downDynamicPseudoRatio_/((double) numberTimesDown_)),
959             shadowEstimateDown_);
960    else
961      printf("Updating id %d - down change %g (true %g) - shadow estimate %g\n",
962             id_,change,data.change_,shadowEstimateDown_);
963#endif
964    numberTimesDown_++;
965    downDynamicPseudoRatio_ += change/shadowEstimateDown_;
966  } else {
967    // up
968    if (!feasible) {
969      double distanceToCutoff=0.0;
970      //double objectiveValue = model_->getCurrentMinimizationObjValue();
971      distanceToCutoff =  model_->getCutoff()  - originalValue;
972      if (distanceToCutoff<1.0e20) 
973        change = distanceToCutoff*2.0;
974      else 
975        change = (upDynamicPseudoRatio_*shadowEstimateUp_+1.0e-3)*10.0;
976    } 
977    change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
978#ifdef PRINT_SHADOW
979    if (numberTimesUp_)
980      printf("Updating id %d - up change %g (true %g) - nup %d estimated change %g - raw shadow estimate %g\n",
981             id_,change,data.change_,numberTimesUp_,shadowEstimateUp_*
982             (upDynamicPseudoRatio_/((double) numberTimesUp_)),
983             shadowEstimateUp_);
984    else
985      printf("Updating id %d - up change %g (true %g) - shadow estimate %g\n",
986             id_,change,data.change_,shadowEstimateUp_);
987#endif
988    numberTimesUp_++;
989    upDynamicPseudoRatio_ += change/shadowEstimateUp_;
990  }
991}
992
993/* Create an OsiSolverBranch object
994   
995This returns NULL if branch not represented by bound changes
996*/
997OsiSolverBranch * 
998CbcSOS::solverBranch() const
999{
1000  int j;
1001  const double * solution = model_->testSolution();
1002  double integerTolerance = 
1003      model_->getDblParam(CbcModel::CbcIntegerTolerance);
1004  OsiSolverInterface * solver = model_->solver();
1005  const double * upper = solver->getColUpper();
1006  int firstNonFixed=-1;
1007  int lastNonFixed=-1;
1008  int firstNonZero=-1;
1009  int lastNonZero = -1;
1010  double weight = 0.0;
1011  double sum =0.0;
1012  double * fix = new double[numberMembers_];
1013  int * which = new int[numberMembers_];
1014  for (j=0;j<numberMembers_;j++) {
1015    int iColumn = members_[j];
1016    // fix all on one side or other (even if fixed)
1017    fix[j]=0.0;
1018    which[j]=iColumn;
1019    if (upper[iColumn]) {
1020      double value = CoinMax(0.0,solution[iColumn]);
1021      sum += value;
1022      if (firstNonFixed<0)
1023        firstNonFixed=j;
1024      lastNonFixed=j;
1025      if (value>integerTolerance) {
1026        weight += weights_[j]*value;
1027        if (firstNonZero<0)
1028          firstNonZero=j;
1029        lastNonZero=j;
1030      }
1031    }
1032  }
1033  assert (lastNonZero-firstNonZero>=sosType_) ;
1034  // find where to branch
1035  assert (sum>0.0);
1036  weight /= sum;
1037  // down branch fixes ones above weight to 0
1038  int iWhere;
1039  int iDownStart=0;
1040  int iUpEnd=0;
1041  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
1042    if (weight<weights_[iWhere+1])
1043      break;
1044  if (sosType_==1) {
1045    // SOS 1
1046    iUpEnd=iWhere+1;
1047    iDownStart=iUpEnd;
1048  } else {
1049    // SOS 2
1050    if (iWhere==firstNonFixed)
1051      iWhere++;;
1052    if (iWhere==lastNonFixed-1)
1053      iWhere = lastNonFixed-2;
1054    iUpEnd=iWhere+1;
1055    iDownStart=iUpEnd+1;
1056  }
1057  //
1058  OsiSolverBranch * branch = new OsiSolverBranch();
1059  branch->addBranch(-1,0,NULL,NULL,numberMembers_-iDownStart,which+iDownStart,fix);
1060  branch->addBranch(1,0,NULL,NULL,iUpEnd,which,fix);
1061  delete [] fix;
1062  delete [] which;
1063  return branch;
1064}
1065// Construct an OsiSOS object
1066OsiSOS * 
1067CbcSOS::osiObject(const OsiSolverInterface * solver) const
1068{
1069  OsiSOS * obj = new OsiSOS(solver,numberMembers_,members_,weights_,sosType_);
1070  obj->setPriority(priority());
1071  return obj;
1072}
1073
1074//##############################################################################
1075
1076/** Default Constructor
1077
1078  Equivalent to an unspecified binary variable.
1079*/
1080CbcSimpleInteger::CbcSimpleInteger ()
1081  : CbcObject(),
1082    originalLower_(0.0),
1083    originalUpper_(1.0),
1084    breakEven_(0.5),
1085    columnNumber_(-1),
1086    preferredWay_(0)
1087{
1088}
1089
1090/** Useful constructor
1091
1092  Loads actual upper & lower bounds for the specified variable.
1093*/
1094CbcSimpleInteger::CbcSimpleInteger ( CbcModel * model, int iColumn, double breakEven)
1095  : CbcObject(model)
1096{
1097  columnNumber_ = iColumn ;
1098  originalLower_ = model->solver()->getColLower()[columnNumber_] ;
1099  originalUpper_ = model->solver()->getColUpper()[columnNumber_] ;
1100  breakEven_ = breakEven;
1101  assert (breakEven_>0.0&&breakEven_<1.0);
1102  preferredWay_ = 0;
1103}
1104
1105
1106// Copy constructor
1107CbcSimpleInteger::CbcSimpleInteger ( const CbcSimpleInteger & rhs)
1108  :CbcObject(rhs)
1109
1110{
1111  columnNumber_ = rhs.columnNumber_;
1112  originalLower_ = rhs.originalLower_;
1113  originalUpper_ = rhs.originalUpper_;
1114  breakEven_ = rhs.breakEven_;
1115  preferredWay_ = rhs.preferredWay_;
1116}
1117
1118// Clone
1119CbcObject *
1120CbcSimpleInteger::clone() const
1121{
1122  return new CbcSimpleInteger(*this);
1123}
1124
1125// Assignment operator
1126CbcSimpleInteger & 
1127CbcSimpleInteger::operator=( const CbcSimpleInteger& rhs)
1128{
1129  if (this!=&rhs) {
1130    CbcObject::operator=(rhs);
1131    columnNumber_ = rhs.columnNumber_;
1132    originalLower_ = rhs.originalLower_;
1133    originalUpper_ = rhs.originalUpper_;
1134    breakEven_ = rhs.breakEven_;
1135    preferredWay_ = rhs.preferredWay_;
1136  }
1137  return *this;
1138}
1139
1140// Destructor
1141CbcSimpleInteger::~CbcSimpleInteger ()
1142{
1143}
1144// Construct an OsiSimpleInteger object
1145OsiSimpleInteger * 
1146CbcSimpleInteger::osiObject() const
1147{
1148  OsiSimpleInteger * obj = new OsiSimpleInteger(columnNumber_,
1149                                                originalLower_,originalUpper_);
1150  obj->setPriority(priority());
1151  return obj;
1152}
1153
1154double
1155CbcSimpleInteger::infeasibility(const OsiSolverInterface * solver, const OsiBranchingInformation * info,
1156                         int & preferredWay) const
1157{
1158  double value = info->solution_[columnNumber_];
1159  value = CoinMax(value, info->lower_[columnNumber_]);
1160  value = CoinMin(value, info->upper_[columnNumber_]);
1161  double nearest = floor(value+(1.0-breakEven_));
1162  assert (breakEven_>0.0&&breakEven_<1.0);
1163  if (nearest>value) 
1164    preferredWay=1;
1165  else
1166    preferredWay=-1;
1167  if (preferredWay_)
1168    preferredWay=preferredWay_;
1169  double weight = fabs(value-nearest);
1170  // normalize so weight is 0.5 at break even
1171  if (nearest<value)
1172    weight = (0.5/breakEven_)*weight;
1173  else
1174    weight = (0.5/(1.0-breakEven_))*weight;
1175  if (fabs(value-nearest)<=info->integerTolerance_) 
1176    return 0.0;
1177  else
1178    return weight;
1179}
1180double 
1181CbcSimpleInteger::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const 
1182{
1183  double value = info->solution_[columnNumber_];
1184#ifdef COIN_DEVELOP
1185  if (fabs(value-floor(value+0.5))>1.0e-5)
1186    printf("value for %d away from integer %g\n",columnNumber_,value);
1187#endif
1188  double newValue = CoinMax(value, info->lower_[columnNumber_]);
1189  newValue = CoinMin(newValue, info->upper_[columnNumber_]);
1190  newValue = floor(newValue+0.5);
1191  solver->setColLower(columnNumber_,newValue);
1192  solver->setColUpper(columnNumber_,newValue);
1193  return fabs(value-newValue);
1194}
1195
1196/* Create an OsiSolverBranch object
1197   
1198This returns NULL if branch not represented by bound changes
1199*/
1200OsiSolverBranch * 
1201CbcSimpleInteger::solverBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
1202{
1203  double value = info->solution_[columnNumber_];
1204  value = CoinMax(value, info->lower_[columnNumber_]);
1205  value = CoinMin(value, info->upper_[columnNumber_]);
1206  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
1207#ifndef NDEBUG
1208  double nearest = floor(value+0.5);
1209  assert (fabs(value-nearest)>info->integerTolerance_);
1210#endif
1211  OsiSolverBranch * branch = new OsiSolverBranch();
1212  branch->addBranch(columnNumber_,value);
1213  return branch;
1214}
1215// Creates a branching object
1216CbcBranchingObject * 
1217CbcSimpleInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) 
1218{
1219  CbcIntegerBranchingObject * branch = new CbcIntegerBranchingObject(model_,0,-1,0.5);
1220  fillCreateBranch(branch,info,way);
1221  return branch;
1222}
1223// Fills in a created branching object
1224void 
1225CbcSimpleInteger::fillCreateBranch(CbcIntegerBranchingObject * branch, const OsiBranchingInformation * info, int way) 
1226{
1227  branch->setOriginalObject(this);
1228  double value = info->solution_[columnNumber_];
1229  value = CoinMax(value, info->lower_[columnNumber_]);
1230  value = CoinMin(value, info->upper_[columnNumber_]);
1231  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
1232  if (!info->hotstartSolution_) {
1233#ifndef NDEBUG
1234    double nearest = floor(value+0.5);
1235    assert (fabs(value-nearest)>info->integerTolerance_);
1236#endif
1237  } else {
1238    double targetValue = info->hotstartSolution_[columnNumber_];
1239    if (way>0)
1240      value = targetValue-0.1;
1241    else
1242      value = targetValue+0.1;
1243  }
1244  branch->fillPart(columnNumber_,way,value);
1245}
1246/* Column number if single column object -1 otherwise,
1247   so returns >= 0
1248   Used by heuristics
1249*/
1250int 
1251CbcSimpleInteger::columnNumber() const
1252{
1253  return columnNumber_;
1254}
1255/* Reset variable bounds to their original values.
1256 
1257    Bounds may be tightened, so it may be good to be able to set this info in object.
1258*/
1259void 
1260CbcSimpleInteger::resetBounds(const OsiSolverInterface * solver) 
1261{
1262  originalLower_ = solver->getColLower()[columnNumber_] ;
1263  originalUpper_ = solver->getColUpper()[columnNumber_] ;
1264}
1265
1266/*  Change column numbers after preprocessing
1267 */
1268void 
1269CbcSimpleInteger::resetSequenceEtc(int numberColumns, const int * originalColumns) 
1270{
1271  assert (numberColumns>0);
1272  int iColumn;
1273#if 0
1274  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1275    if (columnNumber_==originalColumns[iColumn])
1276      break;
1277  }
1278  assert (iColumn<numberColumns);
1279#else
1280  iColumn=originalColumns[columnNumber_];
1281  assert (iColumn>=0);
1282#endif
1283  columnNumber_ = iColumn;
1284}
1285
1286// Infeasibility - large is 0.5
1287double 
1288CbcSimpleInteger::infeasibility(int & preferredWay) const
1289{
1290  OsiBranchingInformation info(model_->solver(),model_->normalSolver(),false);
1291  return infeasibility(model_->solver(),&info,preferredWay);
1292}
1293
1294// This looks at solution and sets bounds to contain solution
1295/** More precisely: it first forces the variable within the existing
1296    bounds, and then tightens the bounds to fix the variable at the
1297    nearest integer value.
1298*/
1299void 
1300CbcSimpleInteger::feasibleRegion()
1301{
1302  abort();
1303}
1304CbcBranchingObject * 
1305CbcSimpleInteger::createBranch( int way) 
1306{
1307  abort();
1308  return NULL;
1309}
1310
1311//##############################################################################
1312
1313// Default Constructor
1314CbcIntegerBranchingObject::CbcIntegerBranchingObject()
1315  :CbcBranchingObject()
1316{
1317  down_[0] = 0.0;
1318  down_[1] = 0.0;
1319  up_[0] = 0.0;
1320  up_[1] = 0.0;
1321#ifdef FUNNY_BRANCHING
1322  variables_ = NULL;
1323  newBounds_ = NULL;
1324  numberExtraChangedBounds_ = 0;
1325#endif
1326}
1327// Useful constructor
1328CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
1329                                                      int variable, int way , double value)
1330  :CbcBranchingObject(model,variable,way,value)
1331{
1332  int iColumn = variable;
1333  assert (model_->solver()->getNumCols()>0);
1334  down_[0] = model_->solver()->getColLower()[iColumn];
1335  down_[1] = floor(value_);
1336  up_[0] = ceil(value_);
1337  up_[1] = model->getColUpper()[iColumn];
1338#ifdef FUNNY_BRANCHING
1339  variables_ = NULL;
1340  newBounds_ = NULL;
1341  numberExtraChangedBounds_ = 0;
1342#endif
1343}
1344// Does part of constructor
1345void 
1346CbcIntegerBranchingObject::fillPart (int variable,
1347                                 int way , double value) 
1348{
1349  //originalObject_=NULL;
1350  branchIndex_=0;
1351  value_=value;
1352  numberBranches_=2;
1353  //model_= model;
1354  //originalCbcObject_=NULL;
1355  variable_=variable;
1356  way_=way;
1357  int iColumn = variable;
1358  down_[0] = model_->solver()->getColLower()[iColumn];
1359  down_[1] = floor(value_);
1360  up_[0] = ceil(value_);
1361  up_[1] = model_->getColUpper()[iColumn];
1362}
1363// Useful constructor for fixing
1364CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
1365                                                      int variable, int way,
1366                                                      double lowerValue, 
1367                                                      double upperValue)
1368  :CbcBranchingObject(model,variable,way,lowerValue)
1369{
1370  setNumberBranchesLeft(1);
1371  down_[0] = lowerValue;
1372  down_[1] = upperValue;
1373  up_[0] = lowerValue;
1374  up_[1] = upperValue;
1375#ifdef FUNNY_BRANCHING
1376  variables_ = NULL;
1377  newBounds_ = NULL;
1378  numberExtraChangedBounds_ = 0;
1379#endif
1380}
1381 
1382
1383// Copy constructor
1384CbcIntegerBranchingObject::CbcIntegerBranchingObject ( const CbcIntegerBranchingObject & rhs) :CbcBranchingObject(rhs)
1385{
1386  down_[0] = rhs.down_[0];
1387  down_[1] = rhs.down_[1];
1388  up_[0] = rhs.up_[0];
1389  up_[1] = rhs.up_[1];
1390#ifdef FUNNY_BRANCHING
1391  numberExtraChangedBounds_ = rhs.numberExtraChangedBounds_;
1392  int size = numberExtraChangedBounds_*(sizeof(double)+sizeof(int));
1393  char * temp = new char [size];
1394  newBounds_ = (double *) temp;
1395  variables_ = (int *) (newBounds_+numberExtraChangedBounds_);
1396
1397  int i ;
1398  for (i=0;i<numberExtraChangedBounds_;i++) {
1399    variables_[i]=rhs.variables_[i];
1400    newBounds_[i]=rhs.newBounds_[i];
1401  }
1402#endif
1403}
1404
1405// Assignment operator
1406CbcIntegerBranchingObject & 
1407CbcIntegerBranchingObject::operator=( const CbcIntegerBranchingObject& rhs)
1408{
1409  if (this != &rhs) {
1410    CbcBranchingObject::operator=(rhs);
1411    down_[0] = rhs.down_[0];
1412    down_[1] = rhs.down_[1];
1413    up_[0] = rhs.up_[0];
1414    up_[1] = rhs.up_[1];
1415#ifdef FUNNY_BRANCHING
1416    delete [] newBounds_;
1417    numberExtraChangedBounds_ = rhs.numberExtraChangedBounds_;
1418    int size = numberExtraChangedBounds_*(sizeof(double)+sizeof(int));
1419    char * temp = new char [size];
1420    newBounds_ = (double *) temp;
1421    variables_ = (int *) (newBounds_+numberExtraChangedBounds_);
1422   
1423    int i ;
1424    for (i=0;i<numberExtraChangedBounds_;i++) {
1425      variables_[i]=rhs.variables_[i];
1426      newBounds_[i]=rhs.newBounds_[i];
1427    }
1428#endif
1429  }
1430  return *this;
1431}
1432CbcBranchingObject * 
1433CbcIntegerBranchingObject::clone() const
1434{ 
1435  return (new CbcIntegerBranchingObject(*this));
1436}
1437
1438
1439// Destructor
1440CbcIntegerBranchingObject::~CbcIntegerBranchingObject ()
1441{
1442  // for debugging threads
1443  way_=-23456789;
1444#ifdef FUNNY_BRANCHING
1445  delete [] newBounds_;
1446#endif
1447}
1448
1449/*
1450  Perform a branch by adjusting the bounds of the specified variable. Note
1451  that each arm of the branch advances the object to the next arm by
1452  advancing the value of way_.
1453
1454  Providing new values for the variable's lower and upper bounds for each
1455  branching direction gives a little bit of additional flexibility and will
1456  be easily extensible to multi-way branching.
1457  Returns change in guessed objective on next branch
1458*/
1459double
1460CbcIntegerBranchingObject::branch()
1461{
1462  // for debugging threads
1463  if (way_<-1||way_>100000) {
1464    printf("way %d, left %d, iCol %d, variable %d\n",
1465           way_,numberBranchesLeft(),
1466           originalCbcObject_->columnNumber(),variable_);
1467    assert (way_!=-23456789);
1468  }
1469  decrementNumberBranchesLeft();
1470  if (down_[1]==-COIN_DBL_MAX)
1471    return 0.0;
1472  int iColumn = originalCbcObject_->columnNumber();
1473  assert (variable_==iColumn);
1474  double olb,oub ;
1475  olb = model_->solver()->getColLower()[iColumn] ;
1476  oub = model_->solver()->getColUpper()[iColumn] ;
1477#ifdef COIN_DEVELOP
1478  if (olb!=down_[0]||oub!=up_[1]) {
1479    if (way_>0)
1480      printf("branching up on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
1481             iColumn,olb,oub,up_[0],up_[1],down_[0],down_[1]) ; 
1482    else
1483      printf("branching down on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
1484             iColumn,olb,oub,down_[0],down_[1],up_[0],up_[1]) ; 
1485  }
1486#endif
1487  if (way_<0) {
1488#ifdef CBC_DEBUG
1489  { double olb,oub ;
1490    olb = model_->solver()->getColLower()[iColumn] ;
1491    oub = model_->solver()->getColUpper()[iColumn] ;
1492    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
1493           iColumn,olb,oub,down_[0],down_[1]) ; }
1494#endif
1495    model_->solver()->setColLower(iColumn,down_[0]);
1496    model_->solver()->setColUpper(iColumn,down_[1]);
1497    //#define CBC_PRINT2
1498#ifdef CBC_PRINT2
1499    printf("%d branching down has bounds %g %g",iColumn,down_[0],down_[1]);
1500#endif
1501#ifdef FUNNY_BRANCHING
1502    // branch - do extra bounds
1503    for (int i=0;i<numberExtraChangedBounds_;i++) {
1504      int variable = variables_[i];
1505      if ((variable&0x40000000)!=0) {
1506        // for going down
1507        int k = variable&0x3fffffff;
1508        assert (k!=iColumn);
1509        if ((variable&0x80000000)==0) {
1510          // lower bound changing
1511#ifdef CBC_PRINT2
1512          printf(" extra for %d changes lower from %g to %g",
1513                 k,model_->solver()->getColLower()[k],newBounds_[i]);
1514#endif
1515          model_->solver()->setColLower(k,newBounds_[i]);
1516        } else {
1517          // upper bound changing
1518#ifdef CBC_PRINT2
1519          printf(" extra for %d changes upper from %g to %g",
1520                 k,model_->solver()->getColUpper()[k],newBounds_[i]);
1521#endif
1522          model_->solver()->setColUpper(k,newBounds_[i]);
1523        }
1524      }
1525    }
1526#endif
1527#ifdef CBC_PRINT2
1528    printf("\n");
1529#endif
1530    way_=1;
1531  } else {
1532#ifdef CBC_DEBUG
1533  { double olb,oub ;
1534    olb = model_->solver()->getColLower()[iColumn] ;
1535    oub = model_->solver()->getColUpper()[iColumn] ;
1536    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
1537           iColumn,olb,oub,up_[0],up_[1]) ; }
1538#endif
1539    model_->solver()->setColLower(iColumn,up_[0]);
1540    model_->solver()->setColUpper(iColumn,up_[1]);
1541#ifdef CBC_PRINT2
1542    printf("%d branching up has bounds %g %g",iColumn,up_[0],up_[1]);
1543#endif
1544#ifdef FUNNY_BRANCHING
1545    // branch - do extra bounds
1546    for (int i=0;i<numberExtraChangedBounds_;i++) {
1547      int variable = variables_[i];
1548      if ((variable&0x40000000)==0) {
1549        // for going up
1550        int k = variable&0x3fffffff;
1551        assert (k!=iColumn);
1552        if ((variable&0x80000000)==0) {
1553          // lower bound changing
1554#ifdef CBC_PRINT2
1555          printf(" extra for %d changes lower from %g to %g",
1556                 k,model_->solver()->getColLower()[k],newBounds_[i]);
1557#endif
1558          model_->solver()->setColLower(k,newBounds_[i]);
1559        } else {
1560          // upper bound changing
1561#ifdef CBC_PRINT2
1562          printf(" extra for %d changes upper from %g to %g",
1563                 k,model_->solver()->getColUpper()[k],newBounds_[i]);
1564#endif
1565          model_->solver()->setColUpper(k,newBounds_[i]);
1566        }
1567      }
1568    }
1569#endif
1570#ifdef CBC_PRINT2
1571    printf("\n");
1572#endif
1573    way_=-1;      // Swap direction
1574  }
1575  double nlb = model_->solver()->getColLower()[iColumn];
1576  double nub = model_->solver()->getColUpper()[iColumn];
1577  if (nlb<olb) {
1578#ifndef NDEBUG
1579    printf("bad lb change for column %d from %g to %g\n",iColumn,olb,nlb);
1580#endif
1581    model_->solver()->setColLower(iColumn,CoinMin(olb,nub));
1582    nlb=olb;
1583  }
1584  if (nub>oub) {
1585#ifndef NDEBUG
1586    printf("bad ub change for column %d from %g to %g\n",iColumn,oub,nub);
1587#endif
1588    model_->solver()->setColUpper(iColumn,CoinMax(oub,nlb));
1589  }
1590#ifndef NDEBUG
1591  if (nlb<olb+1.0e-8&&nub>oub-1.0e-8&&false)
1592    printf("bad null change for column %d - bounds %g,%g\n",iColumn,olb,oub);
1593#endif
1594  return 0.0;
1595}
1596/* Update bounds in solver as in 'branch' and update given bounds.
1597   branchState is -1 for 'down' +1 for 'up' */
1598void 
1599CbcIntegerBranchingObject::fix(OsiSolverInterface * solver,
1600                               double * lower, double * upper,
1601                               int branchState) const 
1602{
1603  int iColumn = originalCbcObject_->columnNumber();
1604  assert (variable_==iColumn);
1605  if (branchState<0) {
1606    model_->solver()->setColLower(iColumn,down_[0]);
1607    lower[iColumn]=down_[0];
1608    model_->solver()->setColUpper(iColumn,down_[1]);
1609    upper[iColumn]=down_[1];
1610  } else {
1611    model_->solver()->setColLower(iColumn,up_[0]);
1612    lower[iColumn]=up_[0];
1613    model_->solver()->setColUpper(iColumn,up_[1]);
1614    upper[iColumn]=up_[1];
1615  }
1616}
1617#ifdef FUNNY_BRANCHING
1618// Deactivate bounds for branching
1619void 
1620CbcIntegerBranchingObject::deactivate()
1621{
1622  down_[1]=-COIN_DBL_MAX;
1623}
1624int
1625CbcIntegerBranchingObject::applyExtraBounds(int iColumn, double lower, double upper, int way)
1626{
1627  // branch - do bounds
1628
1629  int i;
1630  int found=0;
1631  if (variable_==iColumn) {
1632    printf("odd applyExtra %d\n",iColumn);
1633    if (way<0) {
1634      down_[0]=CoinMax(lower,down_[0]);
1635      down_[1]=CoinMin(upper,down_[1]);
1636      assert (down_[0]<=down_[1]);
1637    } else {
1638      up_[0]=CoinMax(lower,up_[0]);
1639      up_[1]=CoinMin(upper,up_[1]);
1640      assert (up_[0]<=up_[1]);
1641    }
1642    return 0;
1643  }
1644  int check = (way<0) ? 0x40000000 : 0;
1645  double newLower=lower;
1646  double newUpper=upper;
1647  for (i=0;i<numberExtraChangedBounds_;i++) {
1648    int variable = variables_[i];
1649    if ((variable&0x40000000)==check) {
1650      int k = variable&0x3fffffff;
1651      if (k==iColumn) {
1652        if ((variable&0x80000000)==0) {
1653          // lower bound changing
1654          found |= 1;
1655          newBounds_[i] = CoinMax(lower,newBounds_[i]);
1656          newLower = newBounds_[i];
1657        } else {
1658          // upper bound changing
1659          found |= 2;
1660          newBounds_[i] = CoinMin(upper,newBounds_[i]);
1661          newUpper = newBounds_[i];
1662        }
1663      }
1664    }
1665  }
1666  int nAdd=0;
1667  if ((found&2)==0) {
1668    // need to add new upper
1669    nAdd++;
1670  }
1671  if ((found&1)==0) {
1672    // need to add new lower
1673    nAdd++;
1674  }
1675  if (nAdd) { 
1676    int size = (numberExtraChangedBounds_+nAdd)*(sizeof(double)+sizeof(int));
1677    char * temp = new char [size];
1678    double * newBounds = (double *) temp;
1679    int * variables = (int *) (newBounds+numberExtraChangedBounds_+nAdd);
1680
1681    int i ;
1682    for (i=0;i<numberExtraChangedBounds_;i++) {
1683      variables[i]=variables_[i];
1684      newBounds[i]=newBounds_[i];
1685    }
1686    delete [] newBounds_;
1687    newBounds_ = newBounds;
1688    variables_ = variables;
1689    if ((found&2)==0) {
1690      // need to add new upper
1691      int variable = iColumn | 0x80000000;
1692      variables_[numberExtraChangedBounds_]=variable;
1693      newBounds_[numberExtraChangedBounds_++]=newUpper;
1694    }
1695    if ((found&1)==0) {
1696      // need to add new lower
1697      int variable = iColumn;
1698      variables_[numberExtraChangedBounds_]=variable;
1699      newBounds_[numberExtraChangedBounds_++]=newLower;
1700    }
1701  }
1702 
1703  return (newUpper>=newLower) ? 0 : 1;
1704}
1705#endif
1706// Print what would happen 
1707void
1708CbcIntegerBranchingObject::print()
1709{
1710  int iColumn = originalCbcObject_->columnNumber();
1711  assert (variable_==iColumn);
1712  if (way_<0) {
1713  { double olb,oub ;
1714    olb = model_->solver()->getColLower()[iColumn] ;
1715    oub = model_->solver()->getColUpper()[iColumn] ;
1716    printf("CbcInteger would branch down on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1717           iColumn,variable_,olb,oub,down_[0],down_[1]) ; }
1718  } else {
1719  { double olb,oub ;
1720    olb = model_->solver()->getColLower()[iColumn] ;
1721    oub = model_->solver()->getColUpper()[iColumn] ;
1722    printf("CbcInteger would branch up on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1723           iColumn,variable_,olb,oub,up_[0],up_[1]) ; }
1724  }
1725}
1726
1727/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1728    same type and must have the same original object, but they may have
1729    different feasible regions.
1730    Return the appropriate CbcRangeCompare value (first argument being the
1731    sub/superset if that's the case). In case of overlap (and if \c
1732    replaceIfOverlap is true) replace the current branching object with one
1733    whose feasible region is the overlap.
1734   */
1735CbcRangeCompare
1736CbcIntegerBranchingObject::compareBranchingObject
1737(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1738{
1739  const CbcIntegerBranchingObject* br =
1740    dynamic_cast<const CbcIntegerBranchingObject*>(brObj);
1741  assert(br);
1742  double* thisBd = way_ < 0 ? down_ : up_;
1743  const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
1744  return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
1745}
1746
1747//##############################################################################
1748
1749/** Default Constructor
1750
1751  Equivalent to an unspecified binary variable.
1752*/
1753CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ()
1754  : CbcSimpleInteger(),
1755    downPseudoCost_(1.0e-5),
1756    upPseudoCost_(1.0e-5),
1757    upDownSeparator_(-1.0),
1758    method_(0)
1759{
1760}
1761
1762/** Useful constructor
1763
1764  Loads actual upper & lower bounds for the specified variable.
1765*/
1766CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1767                                    int iColumn, double breakEven)
1768  : CbcSimpleInteger(model,iColumn,breakEven)
1769{
1770  const double * cost = model->getObjCoefficients();
1771  double costValue = CoinMax(1.0e-5,fabs(cost[iColumn]));
1772  // treat as if will cost what it says up
1773  upPseudoCost_=costValue;
1774  // and balance at breakeven
1775  downPseudoCost_=((1.0-breakEven_)*upPseudoCost_)/breakEven_;
1776  upDownSeparator_ = -1.0;
1777  method_=0;
1778}
1779
1780/** Useful constructor
1781
1782  Loads actual upper & lower bounds for the specified variable.
1783*/
1784CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1785                                    int iColumn, double downPseudoCost,
1786                                                        double upPseudoCost)
1787  : CbcSimpleInteger(model,iColumn)
1788{
1789  downPseudoCost_ = CoinMax(1.0e-10,downPseudoCost);
1790  upPseudoCost_ = CoinMax(1.0e-10,upPseudoCost);
1791  breakEven_ = upPseudoCost_/(upPseudoCost_+downPseudoCost_);
1792  upDownSeparator_ = -1.0;
1793  method_=0;
1794}
1795// Useful constructor - passed and model index and pseudo costs
1796CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model, int dummy,int iColumn, 
1797                                                        double downPseudoCost, double upPseudoCost)
1798{
1799  *this=CbcSimpleIntegerPseudoCost(model,iColumn,downPseudoCost,upPseudoCost);
1800  columnNumber_=iColumn;
1801}
1802
1803// Copy constructor
1804CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ( const CbcSimpleIntegerPseudoCost & rhs)
1805  :CbcSimpleInteger(rhs),
1806   downPseudoCost_(rhs.downPseudoCost_),
1807   upPseudoCost_(rhs.upPseudoCost_),
1808   upDownSeparator_(rhs.upDownSeparator_),
1809   method_(rhs.method_)
1810
1811{
1812}
1813
1814// Clone
1815CbcObject *
1816CbcSimpleIntegerPseudoCost::clone() const
1817{
1818  return new CbcSimpleIntegerPseudoCost(*this);
1819}
1820
1821// Assignment operator
1822CbcSimpleIntegerPseudoCost & 
1823CbcSimpleIntegerPseudoCost::operator=( const CbcSimpleIntegerPseudoCost& rhs)
1824{
1825  if (this!=&rhs) {
1826    CbcSimpleInteger::operator=(rhs);
1827    downPseudoCost_=rhs.downPseudoCost_;
1828    upPseudoCost_=rhs.upPseudoCost_;
1829    upDownSeparator_=rhs.upDownSeparator_;
1830    method_=rhs.method_;
1831  }
1832  return *this;
1833}
1834
1835// Destructor
1836CbcSimpleIntegerPseudoCost::~CbcSimpleIntegerPseudoCost ()
1837{
1838}
1839// Creates a branching object
1840CbcBranchingObject * 
1841CbcSimpleIntegerPseudoCost::createBranch(int way) 
1842{
1843  OsiSolverInterface * solver = model_->solver();
1844  const double * solution = model_->testSolution();
1845  const double * lower = solver->getColLower();
1846  const double * upper = solver->getColUpper();
1847  double value = solution[columnNumber_];
1848  value = CoinMax(value, lower[columnNumber_]);
1849  value = CoinMin(value, upper[columnNumber_]);
1850#ifndef NDEBUG
1851  double nearest = floor(value+0.5);
1852  double integerTolerance = 
1853    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1854  assert (upper[columnNumber_]>lower[columnNumber_]);
1855#endif
1856  if (!model_->hotstartSolution()) {
1857    assert (fabs(value-nearest)>integerTolerance);
1858  } else {
1859    const double * hotstartSolution = model_->hotstartSolution();
1860    double targetValue = hotstartSolution[columnNumber_];
1861    if (way>0)
1862      value = targetValue-0.1;
1863    else
1864      value = targetValue+0.1;
1865  }
1866  CbcIntegerPseudoCostBranchingObject * newObject = 
1867    new CbcIntegerPseudoCostBranchingObject(model_,columnNumber_,way,
1868                                            value);
1869  double up =  upPseudoCost_*(ceil(value)-value);
1870  double down =  downPseudoCost_*(value-floor(value));
1871  double changeInGuessed=up-down;
1872  if (way>0)
1873    changeInGuessed = - changeInGuessed;
1874  changeInGuessed=CoinMax(0.0,changeInGuessed);
1875  //if (way>0)
1876  //changeInGuessed += 1.0e8; // bias to stay up
1877  newObject->setChangeInGuessed(changeInGuessed);
1878  newObject->setOriginalObject(this);
1879  return newObject;
1880}
1881// Infeasibility - large is 0.5
1882double 
1883CbcSimpleIntegerPseudoCost::infeasibility(int & preferredWay) const
1884{
1885  OsiSolverInterface * solver = model_->solver();
1886  const double * solution = model_->testSolution();
1887  const double * lower = solver->getColLower();
1888  const double * upper = solver->getColUpper();
1889  if (upper[columnNumber_]==lower[columnNumber_]) {
1890    // fixed
1891    preferredWay=1;
1892    return 0.0;
1893  }
1894  double value = solution[columnNumber_];
1895  value = CoinMax(value, lower[columnNumber_]);
1896  value = CoinMin(value, upper[columnNumber_]);
1897  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
1898    solution[columnNumber_],upper[columnNumber_]);*/
1899  double nearest = floor(value+0.5);
1900  double integerTolerance = 
1901    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1902  double below = floor(value+integerTolerance);
1903  double above = below+1.0;
1904  if (above>upper[columnNumber_]) {
1905    above=below;
1906    below = above -1;
1907  }
1908  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1909  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1910  if (downCost>=upCost)
1911    preferredWay=1;
1912  else
1913    preferredWay=-1;
1914  // See if up down choice set
1915  if (upDownSeparator_>0.0) {
1916    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
1917  }
1918  if (preferredWay_)
1919    preferredWay=preferredWay_;
1920  if (fabs(value-nearest)<=integerTolerance) {
1921    return 0.0;
1922  } else {
1923    // can't get at model so 1,2 don't make sense
1924    assert(method_<1||method_>2);
1925    if (!method_)
1926      return CoinMin(downCost,upCost);
1927    else
1928      return CoinMax(downCost,upCost);
1929  }
1930}
1931
1932// Return "up" estimate
1933double 
1934CbcSimpleIntegerPseudoCost::upEstimate() const
1935{
1936  OsiSolverInterface * solver = model_->solver();
1937  const double * solution = model_->testSolution();
1938  const double * lower = solver->getColLower();
1939  const double * upper = solver->getColUpper();
1940  double value = solution[columnNumber_];
1941  value = CoinMax(value, lower[columnNumber_]);
1942  value = CoinMin(value, upper[columnNumber_]);
1943  if (upper[columnNumber_]==lower[columnNumber_]) {
1944    // fixed
1945    return 0.0;
1946  }
1947  double integerTolerance = 
1948    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1949  double below = floor(value+integerTolerance);
1950  double above = below+1.0;
1951  if (above>upper[columnNumber_]) {
1952    above=below;
1953    below = above -1;
1954  }
1955  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1956  return upCost;
1957}
1958// Return "down" estimate
1959double 
1960CbcSimpleIntegerPseudoCost::downEstimate() const
1961{
1962  OsiSolverInterface * solver = model_->solver();
1963  const double * solution = model_->testSolution();
1964  const double * lower = solver->getColLower();
1965  const double * upper = solver->getColUpper();
1966  double value = solution[columnNumber_];
1967  value = CoinMax(value, lower[columnNumber_]);
1968  value = CoinMin(value, upper[columnNumber_]);
1969  if (upper[columnNumber_]==lower[columnNumber_]) {
1970    // fixed
1971    return 0.0;
1972  }
1973  double integerTolerance = 
1974    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1975  double below = floor(value+integerTolerance);
1976  double above = below+1.0;
1977  if (above>upper[columnNumber_]) {
1978    above=below;
1979    below = above -1;
1980  }
1981  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1982  return downCost;
1983}
1984
1985//##############################################################################
1986
1987// Default Constructor
1988CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1989  :CbcIntegerBranchingObject()
1990{
1991  changeInGuessed_=1.0e-5;
1992}
1993
1994// Useful constructor
1995CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1996                                                      int variable, int way , double value)
1997  :CbcIntegerBranchingObject(model,variable,way,value)
1998{
1999}
2000// Useful constructor for fixing
2001CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
2002                                                      int variable, int way,
2003                                                      double lowerValue, 
2004                                                      double upperValue)
2005  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
2006{
2007  changeInGuessed_=1.0e100;
2008}
2009 
2010
2011// Copy constructor
2012CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject ( 
2013                                 const CbcIntegerPseudoCostBranchingObject & rhs)
2014  :CbcIntegerBranchingObject(rhs)
2015{
2016  changeInGuessed_ = rhs.changeInGuessed_;
2017}
2018
2019// Assignment operator
2020CbcIntegerPseudoCostBranchingObject & 
2021CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject& rhs)
2022{
2023  if (this != &rhs) {
2024    CbcIntegerBranchingObject::operator=(rhs);
2025    changeInGuessed_ = rhs.changeInGuessed_;
2026  }
2027  return *this;
2028}
2029CbcBranchingObject * 
2030CbcIntegerPseudoCostBranchingObject::clone() const
2031{ 
2032  return (new CbcIntegerPseudoCostBranchingObject(*this));
2033}
2034
2035
2036// Destructor
2037CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
2038{
2039}
2040
2041/*
2042  Perform a branch by adjusting the bounds of the specified variable. Note
2043  that each arm of the branch advances the object to the next arm by
2044  advancing the value of way_.
2045
2046  Providing new values for the variable's lower and upper bounds for each
2047  branching direction gives a little bit of additional flexibility and will
2048  be easily extensible to multi-way branching.
2049  Returns change in guessed objective on next branch
2050*/
2051double
2052CbcIntegerPseudoCostBranchingObject::branch()
2053{
2054  CbcIntegerBranchingObject::branch();
2055  return changeInGuessed_;
2056}
2057
2058/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2059    same type and must have the same original object, but they may have
2060    different feasible regions.
2061    Return the appropriate CbcRangeCompare value (first argument being the
2062    sub/superset if that's the case). In case of overlap (and if \c
2063    replaceIfOverlap is true) replace the current branching object with one
2064    whose feasible region is the overlap.
2065*/
2066CbcRangeCompare
2067CbcIntegerPseudoCostBranchingObject::compareBranchingObject
2068(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2069{
2070  const CbcIntegerPseudoCostBranchingObject* br =
2071    dynamic_cast<const CbcIntegerPseudoCostBranchingObject*>(brObj);
2072  assert(br);
2073  double* thisBd = way_ < 0 ? down_ : up_;
2074  const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
2075  return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
2076}
2077
2078
2079//##############################################################################
2080
2081// Default Constructor
2082CbcCliqueBranchingObject::CbcCliqueBranchingObject()
2083  :CbcBranchingObject()
2084{
2085  clique_ = NULL;
2086  downMask_[0]=0;
2087  downMask_[1]=0;
2088  upMask_[0]=0;
2089  upMask_[1]=0;
2090}
2091
2092// Useful constructor
2093CbcCliqueBranchingObject::CbcCliqueBranchingObject (CbcModel * model,
2094                                                    const CbcClique * clique,
2095                                                    int way ,
2096                                                    int numberOnDownSide, const int * down,
2097                                                    int numberOnUpSide, const int * up)
2098  :CbcBranchingObject(model,clique->id(),way,0.5)
2099{
2100  clique_ = clique;
2101  downMask_[0]=0;
2102  downMask_[1]=0;
2103  upMask_[0]=0;
2104  upMask_[1]=0;
2105  int i;
2106  for (i=0;i<numberOnDownSide;i++) {
2107    int sequence = down[i];
2108    int iWord = sequence>>5;
2109    int iBit = sequence - 32*iWord;
2110    unsigned int k = 1<<iBit;
2111    downMask_[iWord] |= k;
2112  }
2113  for (i=0;i<numberOnUpSide;i++) {
2114    int sequence = up[i];
2115    int iWord = sequence>>5;
2116    int iBit = sequence - 32*iWord;
2117    unsigned int k = 1<<iBit;
2118    upMask_[iWord] |= k;
2119  }
2120}
2121
2122// Copy constructor
2123CbcCliqueBranchingObject::CbcCliqueBranchingObject ( const CbcCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
2124{
2125  clique_=rhs.clique_;
2126  downMask_[0]=rhs.downMask_[0];
2127  downMask_[1]=rhs.downMask_[1];
2128  upMask_[0]=rhs.upMask_[0];
2129  upMask_[1]=rhs.upMask_[1];
2130}
2131
2132// Assignment operator
2133CbcCliqueBranchingObject & 
2134CbcCliqueBranchingObject::operator=( const CbcCliqueBranchingObject& rhs)
2135{
2136  if (this != &rhs) {
2137    CbcBranchingObject::operator=(rhs);
2138    clique_=rhs.clique_;
2139    downMask_[0]=rhs.downMask_[0];
2140    downMask_[1]=rhs.downMask_[1];
2141    upMask_[0]=rhs.upMask_[0];
2142    upMask_[1]=rhs.upMask_[1];
2143  }
2144  return *this;
2145}
2146CbcBranchingObject * 
2147CbcCliqueBranchingObject::clone() const
2148{ 
2149  return (new CbcCliqueBranchingObject(*this));
2150}
2151
2152
2153// Destructor
2154CbcCliqueBranchingObject::~CbcCliqueBranchingObject ()
2155{
2156}
2157double
2158CbcCliqueBranchingObject::branch()
2159{
2160  decrementNumberBranchesLeft();
2161  int iWord;
2162  int numberMembers = clique_->numberMembers();
2163  const int * which = clique_->members();
2164  const int * integerVariables = model_->integerVariable();
2165  int numberWords=(numberMembers+31)>>5;
2166  // *** for way - up means fix all those in down section
2167  if (way_<0) {
2168#ifdef FULL_PRINT
2169    printf("Down Fix ");
2170#endif
2171    for (iWord=0;iWord<numberWords;iWord++) {
2172      int i;
2173      for (i=0;i<32;i++) {
2174        unsigned int k = 1<<i;
2175        if ((upMask_[iWord]&k)!=0) {
2176          int iColumn = which[i+32*iWord];
2177#ifdef FULL_PRINT
2178          printf("%d ",i+32*iWord);
2179#endif
2180          // fix weak way
2181          if (clique_->type(i+32*iWord))
2182            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
2183          else
2184            model_->solver()->setColLower(integerVariables[iColumn],1.0);
2185        }
2186      }
2187    }
2188    way_=1;       // Swap direction
2189  } else {
2190#ifdef FULL_PRINT
2191    printf("Up Fix ");
2192#endif
2193    for (iWord=0;iWord<numberWords;iWord++) {
2194      int i;
2195      for (i=0;i<32;i++) {
2196        unsigned int k = 1<<i;
2197        if ((downMask_[iWord]&k)!=0) {
2198          int iColumn = which[i+32*iWord];
2199#ifdef FULL_PRINT
2200          printf("%d ",i+32*iWord);
2201#endif
2202          // fix weak way
2203          if (clique_->type(i+32*iWord))
2204            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
2205          else
2206            model_->solver()->setColLower(integerVariables[iColumn],1.0);
2207        }
2208      }
2209    }
2210    way_=-1;      // Swap direction
2211  }
2212#ifdef FULL_PRINT
2213  printf("\n");
2214#endif
2215  return 0.0;
2216}
2217// Print what would happen 
2218void
2219CbcCliqueBranchingObject::print()
2220{
2221  int iWord;
2222  int numberMembers = clique_->numberMembers();
2223  const int * which = clique_->members();
2224  const int * integerVariables = model_->integerVariable();
2225  int numberWords=(numberMembers+31)>>5;
2226  // *** for way - up means fix all those in down section
2227  if (way_<0) {
2228    printf("Clique - Down Fix ");
2229    for (iWord=0;iWord<numberWords;iWord++) {
2230      int i;
2231      for (i=0;i<32;i++) {
2232        unsigned int k = 1<<i;
2233        if ((upMask_[iWord]&k)!=0) {
2234          int iColumn = which[i+32*iWord];
2235          printf("%d ",integerVariables[iColumn]);
2236        }
2237      }
2238    }
2239  } else {
2240    printf("Clique - Up Fix ");
2241    for (iWord=0;iWord<numberWords;iWord++) {
2242      int i;
2243      for (i=0;i<32;i++) {
2244        unsigned int k = 1<<i;
2245        if ((downMask_[iWord]&k)!=0) {
2246          int iColumn = which[i+32*iWord];
2247          printf("%d ",integerVariables[iColumn]);
2248        }
2249      }
2250    }
2251  }
2252  printf("\n");
2253}
2254
2255static inline int
2256CbcCompareCliques(const CbcClique* cl0, const CbcClique* cl1)
2257{
2258  if (cl0->cliqueType() < cl1->cliqueType()) {
2259    return -1;
2260  }
2261  if (cl0->cliqueType() > cl1->cliqueType()) {
2262    return 1;
2263  }
2264  if (cl0->numberMembers() != cl1->numberMembers()) {
2265    return cl0->numberMembers() - cl1->numberMembers();
2266  }
2267  if (cl0->numberNonSOSMembers() != cl1->numberNonSOSMembers()) {
2268    return cl0->numberNonSOSMembers() - cl1->numberNonSOSMembers();
2269  }
2270  return memcmp(cl0->members(), cl1->members(),
2271                cl0->numberMembers() * sizeof(int));
2272}
2273
2274/** Compare the original object of \c this with the original object of \c
2275    brObj. Assumes that there is an ordering of the original objects.
2276    This method should be invoked only if \c this and brObj are of the same
2277    type.
2278    Return negative/0/positive depending on whether \c this is
2279    smaller/same/larger than the argument.
2280*/
2281int
2282CbcCliqueBranchingObject::compareOriginalObject
2283(const CbcBranchingObject* brObj) const
2284{
2285  const CbcCliqueBranchingObject* br =
2286    dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
2287  assert(br);
2288  return CbcCompareCliques(clique_, br->clique_);
2289}
2290
2291/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2292    same type and must have the same original object, but they may have
2293    different feasible regions.
2294    Return the appropriate CbcRangeCompare value (first argument being the
2295    sub/superset if that's the case). In case of overlap (and if \c
2296    replaceIfOverlap is true) replace the current branching object with one
2297    whose feasible region is the overlap.
2298*/
2299CbcRangeCompare
2300CbcCliqueBranchingObject::compareBranchingObject
2301(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2302{
2303  const CbcCliqueBranchingObject* br =
2304    dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
2305  assert(br);
2306  unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
2307  const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
2308  const CoinUInt64 cl0 =
2309    (static_cast<CoinUInt64>(thisMask[0]) << 32) | thisMask[1];
2310  const CoinUInt64 cl1 =
2311    (static_cast<CoinUInt64>(otherMask[0]) << 32) | otherMask[1];
2312  if (cl0 == cl1) {
2313    return CbcRangeSame;
2314  }
2315  const CoinUInt64 cl_intersection = (cl0 & cl1);
2316  if (cl_intersection == cl0) {
2317    return CbcRangeSuperset;
2318  }
2319  if (cl_intersection == cl1) {
2320    return CbcRangeSubset;
2321  }
2322  const CoinUInt64 cl_xor = (cl0 ^ cl1);
2323  if (cl_intersection == 0 && cl_xor == 0) {
2324    return CbcRangeDisjoint;
2325  }
2326  const CoinUInt64 cl_union = (cl0 | cl1);
2327  thisMask[0] = static_cast<unsigned int>(cl_union >> 32);
2328  thisMask[1] = static_cast<unsigned int>(cl_union & 0xffffffff);
2329  return CbcRangeOverlap;
2330}
2331
2332//##############################################################################
2333
2334// Default Constructor
2335CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject()
2336  :CbcBranchingObject()
2337{
2338  clique_=NULL;
2339  downMask_=NULL;
2340  upMask_=NULL;
2341}
2342
2343// Useful constructor
2344CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject (CbcModel * model,
2345                                                            const CbcClique * clique, 
2346                                                            int way ,
2347                                                    int numberOnDownSide, const int * down,
2348                                                    int numberOnUpSide, const int * up)
2349  :CbcBranchingObject(model,clique->id(),way,0.5)
2350{
2351  clique_ = clique;
2352  int numberMembers = clique_->numberMembers();
2353  int numberWords=(numberMembers+31)>>5;
2354  downMask_ = new unsigned int [numberWords];
2355  upMask_ = new unsigned int [numberWords];
2356  memset(downMask_,0,numberWords*sizeof(unsigned int));
2357  memset(upMask_,0,numberWords*sizeof(unsigned int));
2358  int i;
2359  for (i=0;i<numberOnDownSide;i++) {
2360    int sequence = down[i];
2361    int iWord = sequence>>5;
2362    int iBit = sequence - 32*iWord;
2363    unsigned int k = 1<<iBit;
2364    downMask_[iWord] |= k;
2365  }
2366  for (i=0;i<numberOnUpSide;i++) {
2367    int sequence = up[i];
2368    int iWord = sequence>>5;
2369    int iBit = sequence - 32*iWord;
2370    unsigned int k = 1<<iBit;
2371    upMask_[iWord] |= k;
2372  }
2373}
2374
2375// Copy constructor
2376CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject ( const CbcLongCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
2377{
2378  clique_=rhs.clique_;
2379  if (rhs.downMask_) {
2380    int numberMembers = clique_->numberMembers();
2381    int numberWords=(numberMembers+31)>>5;
2382    downMask_ = new unsigned int [numberWords];
2383    memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
2384    upMask_ = new unsigned int [numberWords];
2385    memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
2386  } else {
2387    downMask_=NULL;
2388    upMask_=NULL;
2389  }   
2390}
2391
2392// Assignment operator
2393CbcLongCliqueBranchingObject & 
2394CbcLongCliqueBranchingObject::operator=( const CbcLongCliqueBranchingObject& rhs)
2395{
2396  if (this != &rhs) {
2397    CbcBranchingObject::operator=(rhs);
2398    clique_=rhs.clique_;
2399    delete [] downMask_;
2400    delete [] upMask_;
2401    if (rhs.downMask_) {
2402      int numberMembers = clique_->numberMembers();
2403      int numberWords=(numberMembers+31)>>5;
2404      downMask_ = new unsigned int [numberWords];
2405      memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
2406      upMask_ = new unsigned int [numberWords];
2407      memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
2408    } else {
2409      downMask_=NULL;
2410      upMask_=NULL;
2411    }   
2412  }
2413  return *this;
2414}
2415CbcBranchingObject * 
2416CbcLongCliqueBranchingObject::clone() const
2417{ 
2418  return (new CbcLongCliqueBranchingObject(*this));
2419}
2420
2421
2422// Destructor
2423CbcLongCliqueBranchingObject::~CbcLongCliqueBranchingObject ()
2424{
2425  delete [] downMask_;
2426  delete [] upMask_;
2427}
2428double
2429CbcLongCliqueBranchingObject::branch()
2430{
2431  decrementNumberBranchesLeft();
2432  int iWord;
2433  int numberMembers = clique_->numberMembers();
2434  const int * which = clique_->members();
2435  const int * integerVariables = model_->integerVariable();
2436  int numberWords=(numberMembers+31)>>5;
2437  // *** for way - up means fix all those in down section
2438  if (way_<0) {
2439#ifdef FULL_PRINT
2440    printf("Down Fix ");
2441#endif
2442    for (iWord=0;iWord<numberWords;iWord++) {
2443      int i;
2444      for (i=0;i<32;i++) {
2445        unsigned int k = 1<<i;
2446        if ((upMask_[iWord]&k)!=0) {
2447          int iColumn = which[i+32*iWord];
2448#ifdef FULL_PRINT
2449          printf("%d ",i+32*iWord);
2450#endif
2451          // fix weak way
2452          if (clique_->type(i+32*iWord))
2453            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
2454          else
2455            model_->solver()->setColLower(integerVariables[iColumn],1.0);
2456        }
2457      }
2458    }
2459    way_=1;       // Swap direction
2460  } else {
2461#ifdef FULL_PRINT
2462    printf("Up Fix ");
2463#endif
2464    for (iWord=0;iWord<numberWords;iWord++) {
2465      int i;
2466      for (i=0;i<32;i++) {
2467        unsigned int k = 1<<i;
2468        if ((downMask_[iWord]&k)!=0) {
2469          int iColumn = which[i+32*iWord];
2470#ifdef FULL_PRINT
2471          printf("%d ",i+32*iWord);
2472#endif
2473          // fix weak way
2474          if (clique_->type(i+32*iWord))
2475            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
2476          else
2477            model_->solver()->setColLower(integerVariables[iColumn],1.0);
2478        }
2479      }
2480    }
2481    way_=-1;      // Swap direction
2482  }
2483#ifdef FULL_PRINT
2484  printf("\n");
2485#endif
2486  return 0.0;
2487}
2488void
2489CbcLongCliqueBranchingObject::print()
2490{
2491  int iWord;
2492  int numberMembers = clique_->numberMembers();
2493  const int * which = clique_->members();
2494  const int * integerVariables = model_->integerVariable();
2495  int numberWords=(numberMembers+31)>>5;
2496  // *** for way - up means fix all those in down section
2497  if (way_<0) {
2498    printf("Clique - Down Fix ");
2499    for (iWord=0;iWord<numberWords;iWord++) {
2500      int i;
2501      for (i=0;i<32;i++) {
2502        unsigned int k = 1<<i;
2503        if ((upMask_[iWord]&k)!=0) {
2504          int iColumn = which[i+32*iWord];
2505          printf("%d ",integerVariables[iColumn]);
2506        }
2507      }
2508    }
2509  } else {
2510    printf("Clique - Up Fix ");
2511    for (iWord=0;iWord<numberWords;iWord++) {
2512      int i;
2513      for (i=0;i<32;i++) {
2514        unsigned int k = 1<<i;
2515        if ((downMask_[iWord]&k)!=0) {
2516          int iColumn = which[i+32*iWord];
2517          printf("%d ",integerVariables[iColumn]);
2518        }
2519      }
2520    }
2521  }
2522  printf("\n");
2523}
2524
2525/** Compare the original object of \c this with the original object of \c
2526    brObj. Assumes that there is an ordering of the original objects.
2527    This method should be invoked only if \c this and brObj are of the same
2528    type.
2529    Return negative/0/positive depending on whether \c this is
2530    smaller/same/larger than the argument.
2531*/
2532int
2533CbcLongCliqueBranchingObject::compareOriginalObject
2534(const CbcBranchingObject* brObj) const
2535{
2536  const CbcLongCliqueBranchingObject* br =
2537    dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
2538  assert(br);
2539  return CbcCompareCliques(clique_, br->clique_);
2540}
2541
2542/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2543    same type and must have the same original object, but they may have
2544    different feasible regions.
2545    Return the appropriate CbcRangeCompare value (first argument being the
2546    sub/superset if that's the case). In case of overlap (and if \c
2547    replaceIfOverlap is true) replace the current branching object with one
2548    whose feasible region is the overlap.
2549*/
2550CbcRangeCompare
2551CbcLongCliqueBranchingObject::compareBranchingObject
2552(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2553{
2554  const CbcLongCliqueBranchingObject* br =
2555    dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
2556  assert(br);
2557  const int numberMembers = clique_->numberMembers();
2558  const int numberWords=(numberMembers+31)>>5;
2559  unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
2560  const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
2561
2562  if (memcmp(thisMask, otherMask, numberWords * sizeof(unsigned int)) == 0) {
2563    return CbcRangeSame;
2564  }
2565  bool canBeSuperset = true;
2566  bool canBeSubset = true;
2567  int i;
2568  for (i = numberWords-1; i >= 0 && (canBeSuperset || canBeSubset); --i) {
2569    const unsigned int both = (thisMask[i] & otherMask[i]);
2570    canBeSuperset &= (both == thisMask[i]);
2571    canBeSubset &= (both == otherMask[i]);
2572  }
2573  if (canBeSuperset) {
2574    return CbcRangeSuperset;
2575  }
2576  if (canBeSubset) {
2577    return CbcRangeSubset;
2578  }
2579
2580  for (i = numberWords-1; i >= 0; --i) {
2581    if ((thisMask[i] ^ otherMask[i]) != 0) {
2582      break;
2583    }
2584  }
2585  if (i == -1) { // complement
2586    return CbcRangeDisjoint;
2587  }
2588  // must be overlap
2589  for (i = numberWords-1; i >= 0; --i) {
2590    thisMask[i] |= otherMask[i];
2591  }
2592  return CbcRangeOverlap;
2593}
2594
2595//##############################################################################
2596
2597// Default Constructor
2598CbcSOSBranchingObject::CbcSOSBranchingObject()
2599  :CbcBranchingObject(),
2600   firstNonzero_(-1),
2601   lastNonzero_(-1)
2602{
2603  set_ = NULL;
2604  separator_=0.0;
2605}
2606
2607// Useful constructor
2608CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
2609                                              const CbcSOS * set,
2610                                              int way ,
2611                                              double separator)
2612  :CbcBranchingObject(model,set->id(),way,0.5)
2613{
2614  set_ = set;
2615  separator_ = separator;
2616  computeNonzeroRange();
2617}
2618
2619// Copy constructor
2620CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
2621 :CbcBranchingObject(rhs),
2622  firstNonzero_(rhs.firstNonzero_),
2623  lastNonzero_(rhs.lastNonzero_)
2624{
2625  set_=rhs.set_;
2626  separator_ = rhs.separator_;
2627}
2628
2629// Assignment operator
2630CbcSOSBranchingObject & 
2631CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject& rhs)
2632{
2633  if (this != &rhs) {
2634    CbcBranchingObject::operator=(rhs);
2635    set_=rhs.set_;
2636    separator_ = rhs.separator_;
2637    firstNonzero_ = rhs.firstNonzero_;
2638    lastNonzero_ = rhs.lastNonzero_;
2639  }
2640  return *this;
2641}
2642CbcBranchingObject * 
2643CbcSOSBranchingObject::clone() const
2644{ 
2645  return (new CbcSOSBranchingObject(*this));
2646}
2647
2648
2649// Destructor
2650CbcSOSBranchingObject::~CbcSOSBranchingObject ()
2651{
2652}
2653
2654void
2655CbcSOSBranchingObject::computeNonzeroRange()
2656{
2657  const int numberMembers = set_->numberMembers();
2658  const double * weights = set_->weights();
2659  int i = 0;
2660  if (way_ < 0) {
2661    for ( i=0;i<numberMembers;i++) {
2662      if (weights[i] > separator_)
2663        break;
2664    }
2665    assert (i<numberMembers);
2666    firstNonzero_ = 0;
2667    lastNonzero_ = i;
2668  } else {
2669    for ( i=0;i<numberMembers;i++) {
2670      if (weights[i] >= separator_)
2671        break;
2672    }
2673    assert (i<numberMembers);
2674    firstNonzero_ = i;
2675    lastNonzero_ = numberMembers;
2676  }
2677}
2678
2679double
2680CbcSOSBranchingObject::branch()
2681{
2682  decrementNumberBranchesLeft();
2683  int numberMembers = set_->numberMembers();
2684  const int * which = set_->members();
2685  const double * weights = set_->weights();
2686  OsiSolverInterface * solver = model_->solver();
2687  //const double * lower = solver->getColLower();
2688  //const double * upper = solver->getColUpper();
2689  // *** for way - up means fix all those in down section
2690  if (way_<0) {
2691    int i;
2692    for ( i=0;i<numberMembers;i++) {
2693      if (weights[i] > separator_)
2694        break;
2695    }
2696    assert (i<numberMembers);
2697    for (;i<numberMembers;i++) 
2698      solver->setColUpper(which[i],0.0);
2699    way_=1;       // Swap direction
2700  } else {
2701    int i;
2702    for ( i=0;i<numberMembers;i++) {
2703      if (weights[i] >= separator_)
2704        break;
2705      else
2706        solver->setColUpper(which[i],0.0);
2707    }
2708    assert (i<numberMembers);
2709    way_=-1;      // Swap direction
2710  }
2711  computeNonzeroRange();
2712  return 0.0;
2713}
2714/* Update bounds in solver as in 'branch' and update given bounds.
2715   branchState is -1 for 'down' +1 for 'up' */
2716void 
2717CbcSOSBranchingObject::fix(OsiSolverInterface * solver,
2718                               double * lower, double * upper,
2719                               int branchState) const 
2720{
2721  int numberMembers = set_->numberMembers();
2722  const int * which = set_->members();
2723  const double * weights = set_->weights();
2724  //const double * lower = solver->getColLower();
2725  //const double * upper = solver->getColUpper();
2726  // *** for way - up means fix all those in down section
2727  if (branchState<0) {
2728    int i;
2729    for ( i=0;i<numberMembers;i++) {
2730      if (weights[i] > separator_)
2731        break;
2732    }
2733    assert (i<numberMembers);
2734    for (;i<numberMembers;i++) {
2735      solver->setColUpper(which[i],0.0);
2736      upper[which[i]]=0.0;
2737    }
2738  } else {
2739    int i;
2740    for ( i=0;i<numberMembers;i++) {
2741      if (weights[i] >= separator_) {
2742        break;
2743      } else {
2744        solver->setColUpper(which[i],0.0);
2745        upper[which[i]]=0.0;
2746      }
2747    }
2748    assert (i<numberMembers);
2749  }
2750}
2751// Print what would happen 
2752void
2753CbcSOSBranchingObject::print()
2754{
2755  int numberMembers = set_->numberMembers();
2756  const int * which = set_->members();
2757  const double * weights = set_->weights();
2758  OsiSolverInterface * solver = model_->solver();
2759  //const double * lower = solver->getColLower();
2760  const double * upper = solver->getColUpper();
2761  int first=numberMembers;
2762  int last=-1;
2763  int numberFixed=0;
2764  int numberOther=0;
2765  int i;
2766  for ( i=0;i<numberMembers;i++) {
2767    double bound = upper[which[i]];
2768    if (bound) {
2769      first = CoinMin(first,i);
2770      last = CoinMax(last,i);
2771    }
2772  }
2773  // *** for way - up means fix all those in down section
2774  if (way_<0) {
2775    printf("SOS Down");
2776    for ( i=0;i<numberMembers;i++) {
2777      double bound = upper[which[i]];
2778      if (weights[i] > separator_)
2779        break;
2780      else if (bound)
2781        numberOther++;
2782    }
2783    assert (i<numberMembers);
2784    for (;i<numberMembers;i++) {
2785      double bound = upper[which[i]];
2786      if (bound)
2787        numberFixed++;
2788    }
2789  } else {
2790    printf("SOS Up");
2791    for ( i=0;i<numberMembers;i++) {
2792      double bound = upper[which[i]];
2793      if (weights[i] >= separator_)
2794        break;
2795      else if (bound)
2796        numberFixed++;
2797    }
2798    assert (i<numberMembers);
2799    for (;i<numberMembers;i++) {
2800      double bound = upper[which[i]];
2801      if (bound)
2802        numberOther++;
2803    }
2804  }
2805  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
2806         separator_,which[first],weights[first],which[last],weights[last],numberFixed,numberOther);
2807}
2808 
2809/** Compare the original object of \c this with the original object of \c
2810    brObj. Assumes that there is an ordering of the original objects.
2811    This method should be invoked only if \c this and brObj are of the same
2812    type.
2813    Return negative/0/positive depending on whether \c this is
2814    smaller/same/larger than the argument.
2815*/
2816int
2817CbcSOSBranchingObject::compareOriginalObject
2818(const CbcBranchingObject* brObj) const
2819{
2820  const CbcSOSBranchingObject* br =
2821    dynamic_cast<const CbcSOSBranchingObject*>(brObj);
2822  assert(br);
2823  const CbcSOS* s0 = set_;
2824  const CbcSOS* s1 = br->set_;
2825  if (s0->sosType() != s1->sosType()) {
2826    return s0->sosType() - s1->sosType();
2827  }
2828  if (s0->numberMembers() != s1->numberMembers()) {
2829    return s0->numberMembers() - s1->numberMembers();
2830  }
2831  const int memberCmp = memcmp(s0->members(), s1->members(),
2832                               s0->numberMembers() * sizeof(int));
2833  if (memberCmp != 0) {
2834    return memberCmp;
2835  }
2836  return memcmp(s0->weights(), s1->weights(),
2837                s0->numberMembers() * sizeof(double));
2838}
2839
2840/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2841    same type and must have the same original object, but they may have
2842    different feasible regions.
2843    Return the appropriate CbcRangeCompare value (first argument being the
2844    sub/superset if that's the case). In case of overlap (and if \c
2845    replaceIfOverlap is true) replace the current branching object with one
2846    whose feasible region is the overlap.
2847*/
2848CbcRangeCompare
2849CbcSOSBranchingObject::compareBranchingObject
2850(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2851{
2852  const CbcSOSBranchingObject* br =
2853    dynamic_cast<const CbcSOSBranchingObject*>(brObj);
2854  assert(br);
2855  if (firstNonzero_ < br->firstNonzero_) {
2856    if (lastNonzero_ >= br->lastNonzero_) {
2857      return CbcRangeSuperset;
2858    } else if (lastNonzero_ <= br->firstNonzero_) {
2859      return CbcRangeDisjoint;
2860    } else {
2861      // overlap
2862      if (replaceIfOverlap) {
2863        firstNonzero_ = br->firstNonzero_;
2864      }
2865      return CbcRangeOverlap;
2866    }
2867  } else if (firstNonzero_ > br->firstNonzero_) {
2868    if (lastNonzero_ <= br->lastNonzero_) {
2869      return CbcRangeSubset;
2870    } else if (firstNonzero_ >= br->lastNonzero_) {
2871      return CbcRangeDisjoint;
2872    } else {
2873      // overlap
2874      if (replaceIfOverlap) {
2875        lastNonzero_ = br->lastNonzero_;
2876      }
2877      return CbcRangeOverlap;
2878    }
2879  } else {
2880    if (lastNonzero_ == br->lastNonzero_) {
2881      return CbcRangeSame;
2882    }
2883    return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
2884  }
2885  return CbcRangeSame; // fake return
2886}
2887
2888//##############################################################################
2889
2890// Default Constructor
2891CbcBranchDefaultDecision::CbcBranchDefaultDecision()
2892  :CbcBranchDecision()
2893{
2894  bestCriterion_ = 0.0;
2895  bestChangeUp_ = 0.0;
2896  bestNumberUp_ = 0;
2897  bestChangeDown_ = 0.0;
2898  bestObject_ = NULL;
2899  bestNumberDown_ = 0;
2900}
2901
2902// Copy constructor
2903CbcBranchDefaultDecision::CbcBranchDefaultDecision (
2904                                    const CbcBranchDefaultDecision & rhs)
2905  :CbcBranchDecision(rhs)
2906{
2907  bestCriterion_ = rhs.bestCriterion_;
2908  bestChangeUp_ = rhs.bestChangeUp_;
2909  bestNumberUp_ = rhs.bestNumberUp_;
2910  bestChangeDown_ = rhs.bestChangeDown_;
2911  bestNumberDown_ = rhs.bestNumberDown_;
2912  bestObject_ = rhs.bestObject_;
2913  model_ = rhs.model_;
2914}
2915
2916CbcBranchDefaultDecision::~CbcBranchDefaultDecision()
2917{
2918}
2919
2920// Clone
2921CbcBranchDecision * 
2922CbcBranchDefaultDecision::clone() const
2923{
2924  return new CbcBranchDefaultDecision(*this);
2925}
2926
2927// Initialize i.e. before start of choosing at a node
2928void 
2929CbcBranchDefaultDecision::initialize(CbcModel * model)
2930{
2931  bestCriterion_ = 0.0;
2932  bestChangeUp_ = 0.0;
2933  bestNumberUp_ = 0;
2934  bestChangeDown_ = 0.0;
2935  bestNumberDown_ = 0;
2936  bestObject_ = NULL;
2937  model_ = model;
2938}
2939
2940
2941/*
2942  Simple default decision algorithm. Compare based on infeasibility (numInfUp,
2943  numInfDn) until a solution is found by search, then switch to change in
2944  objective (changeUp, changeDn). Note that bestSoFar is remembered in
2945  bestObject_, so the parameter bestSoFar is unused.
2946*/
2947
2948int
2949CbcBranchDefaultDecision::betterBranch(CbcBranchingObject * thisOne,
2950                            CbcBranchingObject * bestSoFar,
2951                            double changeUp, int numInfUp,
2952                            double changeDn, int numInfDn)
2953{
2954  bool beforeSolution = cbcModel()->getSolutionCount()==
2955    cbcModel()->getNumberHeuristicSolutions();;
2956  int betterWay=0;
2957  if (beforeSolution) {
2958    if (!bestObject_) {
2959      bestNumberUp_=COIN_INT_MAX;
2960      bestNumberDown_=COIN_INT_MAX;
2961    }
2962    // before solution - choose smallest number
2963    // could add in depth as well
2964    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
2965    if (numInfUp<numInfDn) {
2966      if (numInfUp<bestNumber) {
2967        betterWay = 1;
2968      } else if (numInfUp==bestNumber) {
2969        if (changeUp<bestCriterion_)
2970          betterWay=1;
2971      }
2972    } else if (numInfUp>numInfDn) {
2973      if (numInfDn<bestNumber) {
2974        betterWay = -1;
2975      } else if (numInfDn==bestNumber) {
2976        if (changeDn<bestCriterion_)
2977          betterWay=-1;
2978      }
2979    } else {
2980      // up and down have same number
2981      bool better=false;
2982      if (numInfUp<bestNumber) {
2983        better=true;
2984      } else if (numInfUp==bestNumber) {
2985        if (min(changeUp,changeDn)<bestCriterion_)
2986          better=true;;
2987      }
2988      if (better) {
2989        // see which way
2990        if (changeUp<=changeDn)
2991          betterWay=1;
2992        else
2993          betterWay=-1;
2994      }
2995    }
2996  } else {
2997    if (!bestObject_) {
2998      bestCriterion_=-1.0;
2999    }
3000    // got a solution
3001    if (changeUp<=changeDn) {
3002      if (changeUp>bestCriterion_)
3003        betterWay=1;
3004    } else {
3005      if (changeDn>bestCriterion_)
3006        betterWay=-1;
3007    }
3008  }
3009  if (betterWay) {
3010    bestCriterion_ = CoinMin(changeUp,changeDn);
3011    bestChangeUp_ = changeUp;
3012    bestNumberUp_ = numInfUp;
3013    bestChangeDown_ = changeDn;
3014    bestNumberDown_ = numInfDn;
3015    bestObject_=thisOne;
3016    // See if user is overriding way
3017    if (thisOne->object()&&thisOne->object()->preferredWay())
3018      betterWay = thisOne->object()->preferredWay();
3019  }
3020  return betterWay;
3021}
3022/* Sets or gets best criterion so far */
3023void 
3024CbcBranchDefaultDecision::setBestCriterion(double value)
3025{ 
3026  bestCriterion_ = value;
3027}
3028double 
3029CbcBranchDefaultDecision::getBestCriterion() const
3030{ 
3031  return bestCriterion_;
3032}
3033
3034/* Compare N branching objects. Return index of best
3035   and sets way of branching in chosen object.
3036   
3037   This routine is used only after strong branching.
3038*/
3039
3040int
3041CbcBranchDefaultDecision::bestBranch (CbcBranchingObject ** objects, int numberObjects,
3042                                   int numberUnsatisfied,
3043                                   double * changeUp, int * numberInfeasibilitiesUp,
3044                                   double * changeDown, int * numberInfeasibilitiesDown,
3045                                   double objectiveValue) 
3046{
3047
3048  int bestWay=0;
3049  int whichObject = -1;
3050  if (numberObjects) {
3051    CbcModel * model = cbcModel();
3052    // at continuous
3053    //double continuousObjective = model->getContinuousObjective();
3054    //int continuousInfeasibilities = model->getContinuousInfeasibilities();
3055   
3056    // average cost to get rid of infeasibility
3057    //double averageCostPerInfeasibility =
3058    //(objectiveValue-continuousObjective)/
3059    //(double) (abs(continuousInfeasibilities-numberUnsatisfied)+1);
3060    /* beforeSolution is :
3061       0 - before any solution
3062       n - n heuristic solutions but no branched one
3063       -1 - branched solution found
3064    */
3065    int numberSolutions = model->getSolutionCount();
3066    double cutoff = model->getCutoff();
3067    int method=0;
3068    int i;
3069    if (numberSolutions) {
3070      int numberHeuristic = model->getNumberHeuristicSolutions();
3071      if (numberHeuristic<numberSolutions) {
3072        method = 1;
3073      } else {
3074        method = 2;
3075        // look further
3076        for ( i = 0 ; i < numberObjects ; i++) {
3077          int numberNext = numberInfeasibilitiesUp[i];
3078         
3079          if (numberNext<numberUnsatisfied) {
3080            int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
3081            double perUnsatisfied = changeUp[i]/(double) numberUp;
3082            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3083            if (estimatedObjective<cutoff) 
3084              method=3;
3085          }
3086          numberNext = numberInfeasibilitiesDown[i];
3087          if (numberNext<numberUnsatisfied) {
3088            int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
3089            double perUnsatisfied = changeDown[i]/(double) numberDown;
3090            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3091            if (estimatedObjective<cutoff) 
3092              method=3;
3093          }
3094        }
3095      }
3096      method=2;
3097    } else {
3098      method = 0;
3099    }
3100    // Uncomment next to force method 4
3101    //method=4;
3102    /* Methods :
3103       0 - fewest infeasibilities
3104       1 - largest min change in objective
3105       2 - as 1 but use sum of changes if min close
3106       3 - predicted best solution
3107       4 - take cheapest up branch if infeasibilities same
3108    */
3109    int bestNumber=COIN_INT_MAX;
3110    double bestCriterion=-1.0e50;
3111    double alternativeCriterion = -1.0;
3112    double bestEstimate = 1.0e100;
3113    switch (method) {
3114    case 0:
3115      // could add in depth as well
3116      for ( i = 0 ; i < numberObjects ; i++) {
3117        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
3118        if (thisNumber<=bestNumber) {
3119          int betterWay=0;
3120          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
3121            if (numberInfeasibilitiesUp[i]<bestNumber) {
3122              betterWay = 1;
3123            } else {
3124              if (changeUp[i]<bestCriterion)
3125                betterWay=1;
3126            }
3127          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
3128            if (numberInfeasibilitiesDown[i]<bestNumber) {
3129              betterWay = -1;
3130            } else {
3131              if (changeDown[i]<bestCriterion)
3132                betterWay=-1;
3133            }
3134          } else {
3135            // up and down have same number
3136            bool better=false;
3137            if (numberInfeasibilitiesUp[i]<bestNumber) {
3138              better=true;
3139            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
3140              if (min(changeUp[i],changeDown[i])<bestCriterion)
3141                better=true;;
3142            }
3143            if (better) {
3144              // see which way
3145              if (changeUp[i]<=changeDown[i])
3146                betterWay=1;
3147              else
3148                betterWay=-1;
3149            }
3150          }
3151          if (betterWay) {
3152            bestCriterion = min(changeUp[i],changeDown[i]);
3153            bestNumber = thisNumber;
3154            whichObject = i;
3155            bestWay = betterWay;
3156          }
3157        }
3158      }
3159      break;
3160    case 1:
3161      for ( i = 0 ; i < numberObjects ; i++) {
3162        int betterWay=0;
3163        if (changeUp[i]<=changeDown[i]) {
3164          if (changeUp[i]>bestCriterion)
3165            betterWay=1;
3166        } else {
3167          if (changeDown[i]>bestCriterion)
3168            betterWay=-1;
3169        }
3170        if (betterWay) {
3171          bestCriterion = min(changeUp[i],changeDown[i]);
3172          whichObject = i;
3173          bestWay = betterWay;
3174        }
3175      }
3176      break;
3177    case 2:
3178      for ( i = 0 ; i < numberObjects ; i++) {
3179        double change = min(changeUp[i],changeDown[i]);
3180        double sum = changeUp[i] + changeDown[i];
3181        bool take=false;
3182        if (change>1.1*bestCriterion) 
3183          take=true;
3184        else if (change>0.9*bestCriterion&&sum+change>bestCriterion+alternativeCriterion) 
3185          take=true;
3186        if (take) {
3187          if (changeUp[i]<=changeDown[i]) {
3188            if (changeUp[i]>bestCriterion)
3189              bestWay=1;
3190          } else {
3191            if (changeDown[i]>bestCriterion)
3192              bestWay=-1;
3193          }
3194          bestCriterion = change;
3195          alternativeCriterion = sum;
3196          whichObject = i;
3197        }
3198      }
3199      break;
3200    case 3:
3201      for ( i = 0 ; i < numberObjects ; i++) {
3202        int numberNext = numberInfeasibilitiesUp[i];
3203       
3204        if (numberNext<numberUnsatisfied) {
3205          int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
3206          double perUnsatisfied = changeUp[i]/(double) numberUp;
3207          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3208          if (estimatedObjective<bestEstimate) {
3209            bestEstimate = estimatedObjective;
3210            bestWay=1;
3211            whichObject=i;
3212          }
3213        }
3214        numberNext = numberInfeasibilitiesDown[i];
3215        if (numberNext<numberUnsatisfied) {
3216          int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
3217          double perUnsatisfied = changeDown[i]/(double) numberDown;
3218          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
3219          if (estimatedObjective<bestEstimate) {
3220            bestEstimate = estimatedObjective;
3221            bestWay=-1;
3222            whichObject=i;
3223          }
3224        }
3225      }
3226      break;
3227    case 4:
3228      // if number infeas same then cheapest up
3229      // first get best number or when going down
3230      // now choose smallest change up amongst equal number infeas
3231      for ( i = 0 ; i < numberObjects ; i++) {
3232        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
3233        if (thisNumber<=bestNumber) {
3234          int betterWay=0;
3235          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
3236            if (numberInfeasibilitiesUp[i]<bestNumber) {
3237              betterWay = 1;
3238            } else {
3239              if (changeUp[i]<bestCriterion)
3240                betterWay=1;
3241            }
3242          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
3243            if (numberInfeasibilitiesDown[i]<bestNumber) {
3244              betterWay = -1;
3245            } else {
3246              if (changeDown[i]<bestCriterion)
3247                betterWay=-1;
3248            }
3249          } else {
3250            // up and down have same number
3251            bool better=false;
3252            if (numberInfeasibilitiesUp[i]<bestNumber) {
3253              better=true;
3254            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
3255              if (min(changeUp[i],changeDown[i])<bestCriterion)
3256                better=true;;
3257            }
3258            if (better) {
3259              // see which way
3260              if (changeUp[i]<=changeDown[i])
3261                betterWay=1;
3262              else
3263                betterWay=-1;
3264            }
3265          }
3266          if (betterWay) {
3267            bestCriterion = min(changeUp[i],changeDown[i]);
3268            bestNumber = thisNumber;
3269            whichObject = i;
3270            bestWay = betterWay;
3271          }
3272        }
3273      }
3274      bestCriterion=1.0e50;
3275      for ( i = 0 ; i < numberObjects ; i++) {
3276        int thisNumber = numberInfeasibilitiesUp[i];
3277        if (thisNumber==bestNumber&&changeUp) {
3278          if (changeUp[i]<bestCriterion) {
3279            bestCriterion = changeUp[i];
3280            whichObject = i;
3281            bestWay = 1;
3282          }
3283        }
3284      }
3285      break;
3286    }
3287    // set way in best
3288    if (whichObject>=0) {
3289      CbcBranchingObject * bestObject = objects[whichObject];
3290      if (bestObject->object()&&bestObject->object()->preferredWay()) 
3291        bestWay = bestObject->object()->preferredWay();
3292      bestObject->way(bestWay);
3293    } else {
3294      printf("debug\n");
3295    }
3296  }
3297  return whichObject;
3298}
3299
3300//##############################################################################
3301
3302// Default Constructor
3303CbcFollowOn::CbcFollowOn ()
3304  : CbcObject(),
3305    rhs_(NULL)
3306{
3307}
3308
3309// Useful constructor
3310CbcFollowOn::CbcFollowOn (CbcModel * model)
3311  : CbcObject(model)
3312{
3313  assert (model);
3314  OsiSolverInterface * solver = model_->solver();
3315  matrix_ = *solver->getMatrixByCol();
3316  matrix_.removeGaps();
3317  matrix_.setExtraGap(0.0);
3318  matrixByRow_ = *solver->getMatrixByRow();
3319  int numberRows = matrix_.getNumRows();
3320 
3321  rhs_ = new int[numberRows];
3322  int i;
3323  const double * rowLower = solver->getRowLower();
3324  const double * rowUpper = solver->getRowUpper();
3325  // Row copy
3326  const double * elementByRow = matrixByRow_.getElements();
3327  const int * column = matrixByRow_.getIndices();
3328  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
3329  const int * rowLength = matrixByRow_.getVectorLengths();
3330  for (i=0;i<numberRows;i++) {
3331    rhs_[i]=0;
3332    double value = rowLower[i];
3333    if (value==rowUpper[i]) {
3334      if (floor(value)==value&&value>=1.0&&value<10.0) {
3335        // check elements
3336        bool good=true;
3337        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
3338          int iColumn = column[j];
3339          if (!solver->isBinary(iColumn))
3340            good=false;
3341          double elValue = elementByRow[j];
3342          if (floor(elValue)!=elValue||value<1.0)
3343            good=false;
3344        }
3345        if (good)
3346          rhs_[i]=(int) value;
3347      }
3348    }
3349  }
3350}
3351
3352// Copy constructor
3353CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
3354  :CbcObject(rhs),
3355   matrix_(rhs.matrix_),
3356   matrixByRow_(rhs.matrixByRow_)
3357{
3358  int numberRows = matrix_.getNumRows();
3359  rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
3360}
3361
3362// Clone
3363CbcObject *
3364CbcFollowOn::clone() const
3365{
3366  return new CbcFollowOn(*this);
3367}
3368
3369// Assignment operator
3370CbcFollowOn & 
3371CbcFollowOn::operator=( const CbcFollowOn& rhs)
3372{
3373  if (this!=&rhs) {
3374    CbcObject::operator=(rhs);
3375    delete [] rhs_;
3376    matrix_ = rhs.matrix_;
3377    matrixByRow_ = rhs.matrixByRow_;
3378    int numberRows = matrix_.getNumRows();
3379    rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
3380  }
3381  return *this;
3382}
3383
3384// Destructor
3385CbcFollowOn::~CbcFollowOn ()
3386{
3387  delete [] rhs_;
3388}
3389// As some computation is needed in more than one place - returns row
3390int 
3391CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
3392{
3393  int whichRow=-1;
3394  otherRow=-1;
3395  int numberRows = matrix_.getNumRows();
3396 
3397  int i;
3398  // For sorting
3399  int * sort = new int [numberRows];
3400  int * isort = new int [numberRows];
3401  // Column copy
3402  //const double * element = matrix_.getElements();
3403  const int * row = matrix_.getIndices();
3404  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
3405  const int * columnLength = matrix_.getVectorLengths();
3406  // Row copy
3407  const double * elementByRow = matrixByRow_.getElements();
3408  const int * column = matrixByRow_.getIndices();
3409  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
3410  const int * rowLength = matrixByRow_.getVectorLengths();
3411  OsiSolverInterface * solver = model_->solver();
3412  const double * columnLower = solver->getColLower();
3413  const double * columnUpper = solver->getColUpper();
3414  const double * solution = solver->getColSolution();
3415  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
3416  int nSort=0;
3417  for (i=0;i<numberRows;i++) {
3418    if (rhs_[i]) {
3419      // check elements
3420      double smallest=1.0e10;
3421      double largest=0.0;
3422      int rhsValue=rhs_[i];
3423      int number1=0;
3424      int numberUnsatisfied=0;
3425      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
3426        int iColumn = column[j];
3427        double value = elementByRow[j];
3428        double solValue = solution[iColumn];
3429        if (columnLower[iColumn]!=columnUpper[iColumn]) {
3430          smallest = CoinMin(smallest,value);
3431          largest = CoinMax(largest,value);
3432          if (value==1.0)
3433            number1++;
3434          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
3435            numberUnsatisfied++;
3436        } else {
3437          rhsValue -= (int)(value*floor(solValue+0.5));
3438        }
3439      }
3440      if (numberUnsatisfied>1) {
3441        if (smallest<largest) {
3442          // probably no good but check a few things
3443          assert (largest<=rhsValue);
3444          if (number1==1&&largest==rhsValue)
3445            printf("could fix\n");
3446        } else if (largest==rhsValue) {
3447          sort[nSort]=i;
3448          isort[nSort++]=-numberUnsatisfied;
3449        }
3450      }
3451    }
3452  }
3453  if (nSort>1) {
3454    CoinSort_2(isort,isort+nSort,sort);
3455    CoinZeroN(isort,numberRows);
3456    double * other = new double[numberRows];
3457    CoinZeroN(other,numberRows);
3458    int * which = new int[numberRows];
3459    //#define COUNT
3460#ifndef COUNT
3461    bool beforeSolution = model_->getSolutionCount()==0;
3462#endif
3463    for (int k=0;k<nSort-1;k++) {
3464      i=sort[k];
3465      int numberUnsatisfied = 0;
3466      int n=0;
3467      int j;
3468      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
3469        int iColumn = column[j];
3470        if (columnLower[iColumn]!=columnUpper[iColumn]) {
3471          double solValue = solution[iColumn]-columnLower[iColumn];
3472          if (solValue<1.0-integerTolerance&&solValue>integerTolerance) {
3473            numberUnsatisfied++;
3474            for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
3475              int iRow = row[jj];
3476              if (rhs_[iRow]) {
3477                other[iRow]+=solValue;
3478                if (isort[iRow]) {
3479                  isort[iRow]++;
3480                } else {
3481                  isort[iRow]=1;
3482                  which[n++]=iRow;
3483                }
3484              }
3485            }
3486          }
3487        }
3488      }
3489      double total=0.0;
3490      // Take out row
3491      double sumThis=other[i];
3492      other[i]=0.0;
3493      assert (numberUnsatisfied==isort[i]);
3494      // find one nearest half if solution, one if before solution
3495      int iBest=-1;
3496      double dtarget=0.5*total;
3497#ifdef COUNT
3498      int target = (numberUnsatisfied+1)>>1;
3499      int best=numberUnsatisfied;
3500#else
3501      double best;
3502      if (beforeSolution)
3503        best=dtarget;
3504      else
3505        best=1.0e30;
3506#endif
3507      for (j=0;j<n;j++) {
3508        int iRow = which[j];
3509        double dvalue=other[iRow];
3510        other[iRow]=0.0;
3511#ifdef COUNT
3512        int value = isort[iRow];
3513#endif
3514        isort[iRow]=0;
3515        if (fabs(dvalue)<1.0e-8||fabs(sumThis-dvalue)<1.0e-8)
3516          continue;
3517        if (dvalue<integerTolerance||dvalue>1.0-integerTolerance)
3518          continue;
3519#ifdef COUNT
3520        if (abs(value-target)<best&&value!=numberUnsatisfied) {
3521          best=abs(value-target);
3522          iBest=iRow;
3523          if (dvalue<dtarget)
3524            preferredWay=1;
3525          else
3526            preferredWay=-1;
3527        }
3528#else
3529        if (beforeSolution) {
3530          if (fabs(dvalue-dtarget)>best) {
3531            best = fabs(dvalue-dtarget);
3532            iBest=iRow;
3533            if (dvalue<dtarget)
3534              preferredWay=1;
3535            else
3536              preferredWay=-1;
3537          }
3538        } else {
3539          if (fabs(dvalue-dtarget)<best) {
3540            best = fabs(dvalue-dtarget);
3541            iBest=iRow;
3542            if (dvalue<dtarget)
3543              preferredWay=1;
3544            else
3545              preferredWay=-1;
3546          }
3547        }
3548#endif
3549      }
3550      if (iBest>=0) {
3551        whichRow=i;
3552        otherRow=iBest;
3553        break;
3554      }
3555    }
3556    delete [] which;
3557    delete [] other;
3558  }
3559  delete [] sort;
3560  delete [] isort;
3561  return whichRow;
3562}
3563
3564// Infeasibility - large is 0.5
3565double 
3566CbcFollowOn::infeasibility(int & preferredWay) const
3567{
3568  int otherRow=0;
3569  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
3570  if (whichRow<0)
3571    return 0.0;
3572  else
3573  return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
3574}
3575
3576// This looks at solution and sets bounds to contain solution
3577void 
3578CbcFollowOn::feasibleRegion()
3579{
3580}
3581
3582
3583// Creates a branching object
3584CbcBranchingObject * 
3585CbcFollowOn::createBranch(int way) 
3586{
3587  int otherRow=0;
3588  int preferredWay;
3589  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
3590  assert(way==preferredWay);
3591  assert (whichRow>=0);
3592  int numberColumns = matrix_.getNumCols();
3593 
3594  // Column copy
3595  //const double * element = matrix_.getElements();
3596  const int * row = matrix_.getIndices();
3597  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
3598  const int * columnLength = matrix_.getVectorLengths();
3599  // Row copy
3600  //const double * elementByRow = matrixByRow_.getElements();
3601  const int * column = matrixByRow_.getIndices();
3602  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
3603  const int * rowLength = matrixByRow_.getVectorLengths();
3604  OsiSolverInterface * solver = model_->solver();
3605  const double * columnLower = solver->getColLower();
3606  const double * columnUpper = solver->getColUpper();
3607  //const double * solution = solver->getColSolution();
3608  int nUp=0;
3609  int nDown=0;
3610  int * upList = new int[numberColumns];
3611  int * downList = new int[numberColumns];
3612  int j;
3613  for (j=rowStart[whichRow];j<rowStart[whichRow]+rowLength[whichRow];j++) {
3614    int iColumn = column[j];
3615    if (columnLower[iColumn]!=columnUpper[iColumn]) {
3616      bool up=true;
3617      for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
3618        int iRow = row[jj];
3619        if (iRow==otherRow) {
3620          up=false;
3621          break;
3622        }
3623      }
3624      if (up)
3625        upList[nUp++]=iColumn;
3626      else
3627        downList[nDown++]=iColumn;
3628    }
3629  }
3630  //printf("way %d\n",way);
3631  // create object
3632  //printf("would fix %d down and %d up\n",nDown,nUp);
3633  CbcBranchingObject * branch
3634     = new CbcFixingBranchingObject(model_,way,
3635                                         nDown,downList,nUp,upList);
3636  delete [] upList;
3637  delete [] downList;
3638  return branch;
3639}
3640
3641//##############################################################################
3642
3643// Default Constructor
3644CbcFixingBranchingObject::CbcFixingBranchingObject()
3645  :CbcBranchingObject()
3646{
3647  numberDown_=0;
3648  numberUp_=0;
3649  downList_=NULL;
3650  upList_=NULL;
3651}
3652
3653// Useful constructor
3654CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
3655                                                    int way ,
3656                                                    int numberOnDownSide, const int * down,
3657                                                    int numberOnUpSide, const int * up)
3658  :CbcBranchingObject(model,0,way,0.5)
3659{
3660  numberDown_=numberOnDownSide;
3661  numberUp_=numberOnUpSide;
3662  downList_ = CoinCopyOfArray(down,numberDown_);
3663  upList_ = CoinCopyOfArray(up,numberUp_);
3664}
3665
3666// Copy constructor
3667CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) :CbcBranchingObject(rhs)
3668{
3669  numberDown_=rhs.numberDown_;
3670  numberUp_=rhs.numberUp_;
3671  downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
3672  upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
3673}
3674
3675// Assignment operator
3676CbcFixingBranchingObject & 
3677CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject& rhs)
3678{
3679  if (this != &rhs) {
3680    CbcBranchingObject::operator=(rhs);
3681    delete [] downList_;
3682    delete [] upList_;
3683    numberDown_=rhs.numberDown_;
3684    numberUp_=rhs.numberUp_;
3685    downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
3686    upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
3687  }
3688  return *this;
3689}
3690CbcBranchingObject * 
3691CbcFixingBranchingObject::clone() const
3692{ 
3693  return (new CbcFixingBranchingObject(*this));
3694}
3695
3696
3697// Destructor
3698CbcFixingBranchingObject::~CbcFixingBranchingObject ()
3699{
3700  delete [] downList_;
3701  delete [] upList_;
3702}
3703double
3704CbcFixingBranchingObject::branch()
3705{
3706  decrementNumberBranchesLeft();
3707  OsiSolverInterface * solver = model_->solver();
3708  const double * columnLower = solver->getColLower();
3709  int i;
3710  // *** for way - up means fix all those in up section
3711  if (way_<0) {
3712#ifdef FULL_PRINT
3713    printf("Down Fix ");
3714#endif
3715    //printf("Down Fix %d\n",numberDown_);
3716    for (i=0;i<numberDown_;i++) {
3717      int iColumn = downList_[i];
3718      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
3719#ifdef FULL_PRINT
3720      printf("Setting bound on %d to lower bound\n",iColumn);
3721#endif
3722    }
3723    way_=1;       // Swap direction
3724  } else {
3725#ifdef FULL_PRINT
3726    printf("Up Fix ");
3727#endif
3728    //printf("Up Fix %d\n",numberUp_);
3729    for (i=0;i<numberUp_;i++) {
3730      int iColumn = upList_[i];
3731      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
3732#ifdef FULL_PRINT
3733      printf("Setting bound on %d to lower bound\n",iColumn);
3734#endif
3735    }
3736    way_=-1;      // Swap direction
3737  }
3738#ifdef FULL_PRINT
3739  printf("\n");
3740#endif
3741  return 0.0;
3742}
3743void
3744CbcFixingBranchingObject::print()
3745{
3746  int i;
3747  // *** for way - up means fix all those in up section
3748  if (way_<0) {
3749    printf("Down Fix ");
3750    for (i=0;i<numberDown_;i++) {
3751      int iColumn = downList_[i];
3752      printf("%d ",iColumn);
3753    }
3754  } else {
3755    printf("Up Fix ");
3756    for (i=0;i<numberUp_;i++) {
3757      int iColumn = upList_[i];
3758      printf("%d ",iColumn);
3759    }
3760  }
3761  printf("\n");
3762}
3763
3764/** Compare the original object of \c this with the original object of \c
3765    brObj. Assumes that there is an ordering of the original objects.
3766    This method should be invoked only if \c this and brObj are of the same
3767    type.
3768    Return negative/0/positive depending on whether \c this is
3769    smaller/same/larger than the argument.
3770*/
3771int
3772CbcFixingBranchingObject::compareOriginalObject
3773(const CbcBranchingObject* brObj) const
3774{
3775  throw("must implement");
3776}
3777
3778/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
3779    same type and must have the same original object, but they may have
3780    different feasible regions.
3781    Return the appropriate CbcRangeCompare value (first argument being the
3782    sub/superset if that's the case). In case of overlap (and if \c
3783    replaceIfOverlap is true) replace the current branching object with one
3784    whose feasible region is the overlap.
3785   */
3786CbcRangeCompare
3787CbcFixingBranchingObject::compareBranchingObject
3788(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
3789{
3790  const CbcFixingBranchingObject* br =
3791    dynamic_cast<const CbcFixingBranchingObject*>(brObj);
3792  assert(br);
3793  // If two FixingBranchingObject's have the same base object then it's pretty
3794  // much guaranteed
3795  throw("must implement");
3796}
3797
3798//##############################################################################
3799
3800// Default Constructor
3801CbcNWay::CbcNWay ()
3802  : CbcObject(),
3803    numberMembers_(0),
3804    members_(NULL),
3805    consequence_(NULL)
3806{
3807}
3808
3809// Useful constructor (which are integer indices)
3810CbcNWay::CbcNWay (CbcModel * model, int numberMembers,
3811                  const int * which, int identifier)
3812  : CbcObject(model)
3813{
3814  id_=identifier;
3815  numberMembers_=numberMembers;
3816  if (numberMembers_) {
3817    members_ = new int[numberMembers_];
3818    memcpy(members_,which,numberMembers_*sizeof(int));
3819  } else {
3820    members_ = NULL;
3821  }
3822  consequence_ = NULL;
3823}
3824
3825// Copy constructor
3826CbcNWay::CbcNWay ( const CbcNWay & rhs)
3827  :CbcObject(rhs)
3828{
3829  numberMembers_ = rhs.numberMembers_;
3830  consequence_ = NULL;
3831  if (numberMembers_) {
3832    members_ = new int[numberMembers_];
3833    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
3834    if (rhs.consequence_) {
3835      consequence_ = new CbcConsequence * [numberMembers_];
3836      for (int i=0;i<numberMembers_;i++) {
3837        if (rhs.consequence_[i])
3838          consequence_[i]= rhs.consequence_[i]->clone();
3839        else
3840          consequence_[i]=NULL;
3841      }
3842    }
3843  } else {
3844    members_ = NULL;
3845  }
3846}
3847
3848// Clone
3849CbcObject *
3850CbcNWay::clone() const
3851{
3852  return new CbcNWay(*this);
3853}
3854
3855// Assignment operator
3856CbcNWay & 
3857CbcNWay::operator=( const CbcNWay& rhs)
3858{
3859  if (this!=&rhs) {
3860    CbcObject::operator=(rhs);
3861    delete [] members_;
3862    numberMembers_ = rhs.numberMembers_;
3863    if (consequence_) {
3864      for (int i=0;i<numberMembers_;i++) 
3865        delete consequence_[i];
3866      delete [] consequence_;
3867      consequence_=NULL;
3868    }
3869    if (numberMembers_) {
3870      members_ = new int[numberMembers_];
3871      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
3872    } else {
3873      members_ = NULL;
3874    }
3875    if (rhs.consequence_) {
3876      consequence_ = new CbcConsequence * [numberMembers_];
3877      for (int i=0;i<numberMembers_;i++) {
3878        if (rhs.consequence_[i])
3879          consequence_[i]= rhs.consequence_[i]->clone();
3880        else
3881          consequence_[i]=NULL;
3882      }
3883    }
3884  }
3885  return *this;
3886}
3887
3888// Destructor
3889CbcNWay::~CbcNWay ()
3890{
3891  delete [] members_;
3892  if (consequence_) {
3893    for (int i=0;i<numberMembers_;i++) 
3894      delete consequence_[i];
3895    delete [] consequence_;
3896  }
3897}
3898// Set up a consequence for a single member
3899void 
3900CbcNWay::setConsequence(int iColumn, const CbcConsequence & consequence)
3901{
3902  if (!consequence_) {
3903    consequence_ = new CbcConsequence * [numberMembers_];
3904    for (int i=0;i<numberMembers_;i++) 
3905      consequence_[i]=NULL;
3906  }
3907  for (int i=0;i<numberMembers_;i++) {
3908    if (members_[i]==iColumn) {
3909      consequence_[i]=consequence.clone();
3910      break;
3911    }
3912  }
3913}
3914
3915// Applies a consequence for a single member
3916void 
3917CbcNWay::applyConsequence(int iSequence, int state) const
3918{
3919  assert (state==-9999||state==9999);
3920  if (consequence_) {
3921    CbcConsequence * consequence = consequence_[iSequence];
3922    if (consequence) 
3923      consequence->applyToSolver(model_->solver(),state);
3924  }
3925}
3926 
3927// Infeasibility - large is 0.5
3928double 
3929CbcNWay::infeasibility(int & preferredWay) const
3930{
3931  int numberUnsatis=0;
3932  int j;
3933  OsiSolverInterface * solver = model_->solver();
3934  const double * solution = model_->testSolution();
3935  const double * lower = solver->getColLower();
3936  const double * upper = solver->getColUpper();
3937  double largestValue=0.0;
3938 
3939  double integerTolerance = 
3940    model_->getDblParam(CbcModel::CbcIntegerTolerance);
3941
3942  for (j=0;j<numberMembers_;j++) {
3943    int iColumn = members_[j];
3944    double value = solution[iColumn];
3945    value = CoinMax(value, lower[iColumn]);
3946    value = CoinMin(value, upper[iColumn]);
3947    double distance = CoinMin(value-lower[iColumn],upper[iColumn]-value);
3948    if (distance>integerTolerance) {
3949      numberUnsatis++;
3950      largestValue = CoinMax(distance,largestValue);
3951    }
3952  }
3953  preferredWay=1;
3954  if (numberUnsatis) {
3955    return largestValue;
3956  } else {
3957    return 0.0; // satisfied
3958  }
3959}
3960
3961// This looks at solution and sets bounds to contain solution
3962void 
3963CbcNWay::feasibleRegion()
3964{
3965  int j;
3966  OsiSolverInterface * solver = model_->solver();
3967  const double * solution = model_->testSolution();
3968  const double * lower = solver->getColLower();
3969  const double * upper = solver->getColUpper();
3970  double integerTolerance = 
3971    model_->getDblParam(CbcModel::CbcIntegerTolerance);
3972  for (j=0;j<numberMembers_;j++) {
3973    int iColumn = members_[j];
3974    double value = solution[iColumn];
3975    value = CoinMax(value, lower[iColumn]);
3976    value = CoinMin(value, upper[iColumn]);
3977    if (value>=upper[iColumn]-integerTolerance) {
3978      solver->setColLower(iColumn,upper[iColumn]);
3979    } else {
3980      assert (value<=lower[iColumn]+integerTolerance);
3981      solver->setColUpper(iColumn,lower[iColumn]);
3982    }
3983  }
3984}
3985// Redoes data when sequence numbers change
3986void 
3987CbcNWay::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
3988{
3989  model_=model;
3990  int n2=0;
3991  for (int j=0;j<numberMembers_;j++) {
3992    int iColumn = members_[j];
3993    int i;
3994    for (i=0;i<numberColumns;i++) {
3995      if (originalColumns[i]==iColumn)
3996        break;
3997    }
3998    if (i<numberColumns) {
3999      members_[n2]=i;
4000      consequence_[n2++]=consequence_[j];
4001    } else {
4002      delete consequence_[j];
4003    }
4004  }
4005  if (n2<numberMembers_) {
4006    printf("** NWay number of members reduced from %d to %d!\n",numberMembers_,n2);
4007    numberMembers_=n2;
4008  }
4009}
4010
4011
4012// Creates a branching object
4013CbcBranchingObject * 
4014CbcNWay::createBranch(int way) 
4015{
4016  int numberFree=0;
4017  int j;
4018
4019  OsiSolverInterface * solver = model_->solver();
4020  const double * solution = model_->testSolution();
4021  const double * lower = solver->getColLower();
4022  const double * upper = solver->getColUpper();
4023  int * list = new int[numberMembers_];
4024  double * sort = new double[numberMembers_];
4025
4026  for (j=0;j<numberMembers_;j++) {
4027    int iColumn = members_[j];
4028    double value = solution[iColumn];
4029    value = CoinMax(value, lower[iColumn]);
4030    value = CoinMin(value, upper[iColumn]);
4031    if (upper[iColumn]>lower[iColumn]) {
4032      double distance = upper[iColumn]-value;
4033      list[numberFree]=j;
4034      sort[numberFree++]=distance;
4035    }
4036  }
4037  assert (numberFree);
4038  // sort
4039  CoinSort_2(sort,sort+numberFree,list);
4040  // create object
4041  CbcBranchingObject * branch;
4042  branch = new CbcNWayBranchingObject(model_,this,numberFree,list);
4043  branch->setOriginalObject(this);
4044  delete [] list;
4045  delete [] sort;
4046  return branch;
4047}
4048 
4049//##############################################################################
4050
4051// Default Constructor
4052CbcNWayBranchingObject::CbcNWayBranchingObject()
4053  :CbcBranchingObject()
4054{
4055  order_=NULL;
4056  object_=NULL;
4057  numberInSet_=0;
4058  way_=0;
4059}
4060
4061// Useful constructor
4062CbcNWayBranchingObject::CbcNWayBranchingObject (CbcModel * model,
4063                                                const CbcNWay * nway, 
4064                                                int number, const int * order)
4065  :CbcBranchingObject(model,nway->id(),-1,0.5)
4066{
4067  numberBranches_ = number;
4068  order_ = new int [number];
4069  object_=nway;
4070  numberInSet_=number;
4071  memcpy(order_,order,number*sizeof(int));
4072}
4073
4074// Copy constructor
4075CbcNWayBranchingObject::CbcNWayBranchingObject ( const CbcNWayBranchingObject & rhs) :CbcBranchingObject(rhs)
4076{
4077  numberInSet_=rhs.numberInSet_;
4078  object_=rhs.object_;
4079  if (numberInSet_) {
4080    order_ = new int [numberInSet_];
4081    memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
4082  } else {
4083    order_=NULL;
4084  }   
4085}
4086
4087// Assignment operator
4088CbcNWayBranchingObject & 
4089CbcNWayBranchingObject::operator=( const CbcNWayBranchingObject& rhs)
4090{
4091  if (this != &rhs) {
4092    CbcBranchingObject::operator=(rhs);
4093    object_=rhs.object_;
4094    delete [] order_;
4095    numberInSet_=rhs.numberInSet_;
4096    if (numberInSet_) {
4097      order_ = new int [numberInSet_];
4098      memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
4099    } else {
4100      order_=NULL;
4101    }   
4102  }
4103  return *this;
4104}
4105CbcBranchingObject * 
4106CbcNWayBranchingObject::clone() const
4107{ 
4108  return (new CbcNWayBranchingObject(*this));
4109}
4110
4111
4112// Destructor
4113CbcNWayBranchingObject::~CbcNWayBranchingObject ()
4114{
4115  delete [] order_;
4116}
4117double
4118CbcNWayBranchingObject::branch()
4119{
4120  int which = branchIndex_;
4121  branchIndex_++;
4122  assert (numberBranchesLeft()>=0);
4123  if (which==0) {
4124    // first branch so way_ may mean something
4125    assert (way_==-1||way_==1);
4126    if (way_==-1)
4127      which++;
4128  } else if (which==1) {
4129    // second branch so way_ may mean something
4130    assert (way_==-1||way_==1);
4131    if (way_==-1)
4132      which--;
4133    // switch way off
4134    way_=0;
4135  }
4136  const double * lower = model_->solver()->getColLower();
4137  const double * upper = model_->solver()->getColUpper();
4138  const int * members = object_->members();
4139  for (int j=0;j<numberInSet_;j++) {
4140    int iSequence = order_[j];
4141    int iColumn = members[iSequence];
4142    if (j!=which) {
4143      model_->solver()->setColUpper(iColumn,lower[iColumn]);
4144      //model_->solver()->setColLower(iColumn,lower[iColumn]);
4145      assert (lower[iColumn]>-1.0e20);
4146      // apply any consequences
4147      object_->applyConsequence(iSequence,-9999);
4148    } else {
4149      model_->solver()->setColLower(iColumn,upper[iColumn]);
4150      //model_->solver()->setColUpper(iColumn,upper[iColumn]);
4151#ifdef FULL_PRINT
4152      printf("Up Fix %d to %g\n",iColumn,upper[iColumn]);
4153#endif
4154      assert (upper[iColumn]<1.0e20);
4155      // apply any consequences
4156      object_->applyConsequence(iSequence,9999);
4157    }
4158  }
4159  return 0.0;
4160}
4161void
4162CbcNWayBranchingObject::print()
4163{
4164  printf("NWay - Up Fix ");
4165  const int * members = object_->members();
4166  for (int j=0;j<way_;j++) {
4167    int iColumn = members[order_[j]];
4168    printf("%d ",iColumn);
4169  }
4170  printf("\n");
4171}
4172
4173/** Compare the original object of \c this with the original object of \c
4174    brObj. Assumes that there is an ordering of the original objects.
4175    This method should be invoked only if \c this and brObj are of the same
4176    type.
4177    Return negative/0/positive depending on whether \c this is
4178    smaller/same/larger than the argument.
4179*/
4180int
4181CbcNWayBranchingObject::compareOriginalObject
4182(const CbcBranchingObject* brObj) const
4183{
4184  throw("must implement");
4185}
4186
4187/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
4188    same type and must have the same original object, but they may have
4189    different feasible regions.
4190    Return the appropriate CbcRangeCompare value (first argument being the
4191    sub/superset if that's the case). In case of overlap (and if \c
4192    replaceIfOverlap is true) replace the current branching object with one
4193    whose feasible region is the overlap.
4194*/
4195CbcRangeCompare
4196CbcNWayBranchingObject::compareBranchingObject
4197(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
4198{
4199  throw("must implement");
4200}
4201
4202//##############################################################################
4203
4204// Default Constructor
4205CbcFixVariable::CbcFixVariable ()
4206  : CbcConsequence(),
4207    numberStates_(0),
4208    states_(NULL),
4209    startLower_(NULL),
4210    startUpper_(NULL),
4211    newBound_(NULL),
4212    variable_(NULL)
4213{
4214}
4215
4216// One useful Constructor
4217CbcFixVariable::CbcFixVariable (int numberStates,const int * states, const int * numberNewLower, 
4218                                const int ** newLowerValue,
4219                                const int ** lowerColumn,
4220                                const int * numberNewUpper, const int ** newUpperValue,
4221                                const int ** upperColumn)
4222  : CbcConsequence(),
4223    states_(NULL),
4224    startLower_(NULL),
4225    startUpper_(NULL),
4226    newBound_(NULL),
4227    variable_(NULL)
4228{
4229  // How much space
4230  numberStates_ = numberStates;
4231  if (numberStates_) {
4232    states_ = new int[numberStates_];
4233    memcpy(states_,states,numberStates_*sizeof(int));
4234    int i;
4235    int n=0;
4236    startLower_ = new int[numberStates_+1];
4237    startUpper_ = new int[numberStates_+1];
4238    startLower_[0]=0;
4239    //count
4240    for (i=0;i<numberStates_;i++) {
4241      n += numberNewLower[i];
4242      startUpper_[i]=n;
4243      n += numberNewUpper[i];
4244      startLower_[i+1]=n;
4245    }
4246    newBound_ = new double [n];
4247    variable_ = new int [n];
4248    n=0;
4249    for (i=0;i<numberStates_;i++) {
4250      int j;
4251      int k;
4252      const int * bound;
4253      const int * variable;
4254      k=numberNewLower[i];
4255      bound = newLowerValue[i];
4256      variable = lowerColumn[i];
4257      for (j=0;j<k;j++) {
4258        newBound_[n]=bound[j];
4259        variable_[n++]=variable[j];
4260      }
4261      k=numberNewUpper[i];
4262      bound = newUpperValue[i];
4263      variable = upperColumn[i];
4264      for (j=0;j<k;j++) {
4265        newBound_[n]=bound[j];
4266        variable_[n++]=variable[j];
4267      }
4268    }
4269  }
4270}
4271
4272// Copy constructor
4273CbcFixVariable::CbcFixVariable ( const CbcFixVariable & rhs)
4274  :CbcConsequence(rhs)
4275{
4276  numberStates_ = rhs.numberStates_;
4277  states_ = NULL;
4278  startLower_ = NULL;
4279  startUpper_ = NULL;
4280  newBound_ = NULL;
4281  variable_ = NULL;
4282  if (numberStates_) {
4283    states_ = CoinCopyOfArray(rhs.states_,numberStates_);
4284    startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
4285    startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
4286    int n=startLower_[numberStates_];
4287    newBound_ = CoinCopyOfArray(rhs.newBound_,n);
4288    variable_ = CoinCopyOfArray(rhs.variable_,n);
4289  }
4290}
4291
4292// Clone
4293CbcConsequence *
4294CbcFixVariable::clone() const
4295{
4296  return new CbcFixVariable(*this);
4297}
4298
4299// Assignment operator
4300CbcFixVariable & 
4301CbcFixVariable::operator=( const CbcFixVariable& rhs)
4302{
4303  if (this!=&rhs) {
4304    CbcConsequence::operator=(rhs);
4305    delete [] states_;
4306    delete [] startLower_;
4307    delete [] startUpper_;
4308    delete [] newBound_;
4309    delete [] variable_;
4310    states_ = NULL;
4311    startLower_ = NULL;
4312    startUpper_ = NULL;
4313    newBound_ = NULL;
4314    variable_ = NULL;
4315    numberStates_ = rhs.numberStates_;
4316    if (numberStates_) {
4317      states_ = CoinCopyOfArray(rhs.states_,numberStates_);
4318      startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
4319      startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
4320      int n=startLower_[numberStates_];
4321      newBound_ = CoinCopyOfArray(rhs.newBound_,n);
4322      variable_ = CoinCopyOfArray(rhs.variable_,n);
4323    }
4324  }
4325  return *this;
4326}
4327
4328// Destructor
4329CbcFixVariable::~CbcFixVariable ()
4330{
4331  delete [] states_;
4332  delete [] startLower_;
4333  delete [] startUpper_;
4334  delete [] newBound_;
4335  delete [] variable_;
4336}
4337// Set up a startLower for a single member
4338void 
4339CbcFixVariable::applyToSolver(OsiSolverInterface * solver, int state) const
4340{
4341  assert (state==-9999||state==9999);
4342  // Find state
4343  int find;
4344  for (find=0;find<numberStates_;find++) 
4345    if (states_[find]==state)
4346      break;
4347  if (find==numberStates_)
4348    return;
4349  int i;
4350  // Set new lower bounds
4351  for (i=startLower_[find];i<startUpper_[find];i++) {
4352    int iColumn = variable_[i];
4353    double value = newBound_[i];
4354    double oldValue = solver->getColLower()[iColumn];
4355    //printf("for %d old lower bound %g, new %g",iColumn,oldValue,value);
4356    solver->setColLower(iColumn,CoinMax(value,oldValue));
4357    //printf(" => %g\n",solver->getColLower()[iColumn]);
4358  }
4359  // Set new upper bounds
4360  for (i=startUpper_[find];i<startLower_[find+1];i++) {
4361    int iColumn = variable_[i];
4362    double value = newBound_[i];
4363    double oldValue = solver->getColUpper()[iColumn];
4364    //printf("for %d old upper bound %g, new %g",iColumn,oldValue,value);
4365    solver->setColUpper(iColumn,CoinMin(value,oldValue));
4366    //printf(" => %g\n",solver->getColUpper()[iColumn]);
4367  }
4368}
4369
4370//##############################################################################
4371
4372// Default Constructor
4373CbcDummyBranchingObject::CbcDummyBranchingObject(CbcModel * model)
4374  :CbcBranchingObject(model,0,0,0.5)
4375{
4376  setNumberBranchesLeft(1);
4377}
4378
4379
4380// Copy constructor
4381CbcDummyBranchingObject::CbcDummyBranchingObject ( const CbcDummyBranchingObject & rhs) :CbcBranchingObject(rhs)
4382{
4383}
4384
4385// Assignment operator
4386CbcDummyBranchingObject & 
4387CbcDummyBranchingObject::operator=( const CbcDummyBranchingObject& rhs)
4388{
4389  if (this != &rhs) {
4390    CbcBranchingObject::operator=(rhs);
4391  }
4392  return *this;
4393}
4394CbcBranchingObject * 
4395CbcDummyBranchingObject::clone() const
4396{ 
4397  return (new CbcDummyBranchingObject(*this));
4398}
4399
4400
4401// Destructor
4402CbcDummyBranchingObject::~CbcDummyBranchingObject ()
4403{
4404}
4405
4406/*
4407  Perform a dummy branch
4408*/
4409double
4410CbcDummyBranchingObject::branch()
4411{
4412  decrementNumberBranchesLeft();
4413  return 0.0;
4414}
4415// Print what would happen 
4416void
4417CbcDummyBranchingObject::print()
4418{
4419  printf("Dummy branch\n");
4420}
4421
4422// Default Constructor
4423CbcGeneral::CbcGeneral() 
4424  : CbcObject()
4425{
4426}
4427
4428// Constructor from model
4429CbcGeneral::CbcGeneral(CbcModel * model)
4430  : CbcObject(model)
4431{
4432}
4433
4434
4435// Destructor
4436CbcGeneral::~CbcGeneral ()
4437{
4438}
4439
4440// Copy constructor
4441CbcGeneral::CbcGeneral ( const CbcGeneral & rhs)
4442  : CbcObject(rhs)
4443{
4444}
4445#ifdef COIN_HAS_CLP
4446#include "OsiClpSolverInterface.hpp"
4447#include "CoinWarmStartBasis.hpp"
4448#include "ClpNode.hpp"
4449#include "CbcBranchDynamic.hpp"
4450// Assignment operator
4451CbcGeneral & 
4452CbcGeneral::operator=( const CbcGeneral& rhs)
4453{
4454  if (this!=&rhs) {
4455    CbcObject::operator=(rhs);
4456  }
4457  return *this;
4458}
4459
4460// Default Constructor
4461CbcGeneralDepth::CbcGeneralDepth ()
4462  : CbcGeneral(),
4463    maximumDepth_(-1),
4464    whichSolution_(-1),
4465    numberNodes_(0),
4466    nodeInfo_(NULL)
4467{
4468}
4469
4470// Useful constructor (which are integer indices)
4471CbcGeneralDepth::CbcGeneralDepth (CbcModel * model, int maximumDepth)
4472  : CbcGeneral(model),
4473    maximumDepth_(maximumDepth),
4474    whichSolution_(-1),
4475    numberNodes_(0),
4476    nodeInfo_(NULL)
4477{
4478  if (maximumDepth_>=0) {
4479    nodeInfo_ = new ClpNodeStuff();
4480    ClpNodeStuff * info = nodeInfo_;
4481    info->nDepth_=maximumDepth_;
4482    int n = (1<<maximumDepth_)+1+maximumDepth_;
4483    ClpNode ** nodeInfo = new ClpNode * [n];
4484    for (int i=0;i<n;i++) 
4485      nodeInfo[i]=NULL;
4486    info->nodeInfo_ = nodeInfo;
4487  } else {
4488    nodeInfo_ = NULL;
4489  }
4490}
4491
4492// Copy constructor
4493CbcGeneralDepth::CbcGeneralDepth ( const CbcGeneralDepth & rhs)
4494  :CbcGeneral(rhs)
4495{
4496  maximumDepth_ = rhs.maximumDepth_;
4497  whichSolution_ = -1;
4498  numberNodes_ = 0;
4499  if (maximumDepth_>=0) {
4500    assert (rhs.nodeInfo_);
4501    nodeInfo_ = new ClpNodeStuff(*rhs.nodeInfo_);
4502    ClpNodeStuff * info = nodeInfo_;
4503    info->nDepth_=maximumDepth_;
4504    if (!info->nodeInfo_) {
4505      int n = (1<<maximumDepth_)+1+maximumDepth_;
4506      ClpNode ** nodeInfo = new ClpNode * [n];
4507      for (int i=0;i<n;i++) 
4508        nodeInfo[i]=NULL;
4509      info->nodeInfo_ = nodeInfo;
4510    }
4511  } else {
4512    nodeInfo_ = NULL;
4513  }
4514}
4515
4516// Clone
4517CbcObject *
4518CbcGeneralDepth::clone() const
4519{
4520  return new CbcGeneralDepth(*this);
4521}
4522
4523// Assignment operator
4524CbcGeneralDepth & 
4525CbcGeneralDepth::operator=( const CbcGeneralDepth& rhs)
4526{
4527  if (this!=&rhs) {
4528    CbcGeneral::operator=(rhs);
4529    delete nodeInfo_;
4530    maximumDepth_ = rhs.maximumDepth_;
4531    whichSolution_ = -1;
4532    numberNodes_ = 0;
4533    if (maximumDepth_>=0) {
4534      assert (rhs.nodeInfo_);
4535      nodeInfo_ = new ClpNodeStuff(*rhs.nodeInfo_);
4536    } else {
4537      nodeInfo_ = NULL;
4538    }
4539  }
4540  return *this;
4541}
4542
4543// Destructor
4544CbcGeneralDepth::~CbcGeneralDepth ()
4545{
4546  delete nodeInfo_;
4547}
4548
4549// Infeasibility - large is 0.5
4550double 
4551CbcGeneralDepth::infeasibility(int & preferredWay) const
4552{
4553  whichSolution_ = -1;
4554  // should use genuine OsiBranchingInformation usefulInfo = model_->usefulInformation();
4555  // for now assume only called when correct
4556  //if (usefulInfo.depth_>=4&&!model_->parentModel()
4557  //     &&(usefulInfo.depth_%2)==0) {
4558  if (true) {
4559    OsiSolverInterface * solver = model_->solver();
4560    OsiClpSolverInterface * clpSolver
4561      = dynamic_cast<OsiClpSolverInterface *> (solver);
4562    if (clpSolver) {
4563      ClpNodeStuff * info = nodeInfo_;
4564      info->integerTolerance_=model_->getIntegerTolerance();
4565      info->integerIncrement_=model_->getCutoffIncrement();
4566      int numberIntegers = model_->numberIntegers();
4567      double * down = new double[numberIntegers];
4568      double * up = new double[numberIntegers];
4569      int * numberDown = new int[numberIntegers];
4570      int * numberUp = new int[numberIntegers];
4571      int * numberDownInfeasible = new int[numberIntegers];
4572      int * numberUpInfeasible = new int[numberIntegers];
4573      model_->fillPseudoCosts(down,up,numberDown,numberUp,
4574                      numberDownInfeasible,numberUpInfeasible);
4575      info->fillPseudoCosts(down,up,numberDown,numberUp,
4576                            numberDownInfeasible,
4577                            numberUpInfeasible,numberIntegers);
4578      info->presolveType_= 1; //1 later
4579      delete [] down;
4580      delete [] up;
4581      delete [] numberDown;
4582      delete [] numberUp;
4583      delete [] numberDownInfeasible;
4584      delete [] numberUpInfeasible;
4585      bool takeHint;
4586      OsiHintStrength strength;
4587      solver->getHintParam(OsiDoReducePrint,takeHint,strength);
4588      ClpSimplex * simplex = clpSolver->getModelPtr();
4589      int saveLevel = simplex->logLevel();
4590      if (strength!=OsiHintIgnore&&takeHint&&saveLevel==1)
4591        simplex->setLogLevel(0);
4592      clpSolver->setBasis();
4593      whichSolution_ = simplex->fathomMany(info);
4594      model_->incrementExtra(info->numberNodesExplored_,
4595                             info->numberIterations_);
4596      // update pseudo costs
4597      double smallest=1.0e50;
4598      double largest=-1.0;
4599      OsiObject ** objects = model_->objects();
4600#ifndef NDEBUG
4601      const int * integerVariable = model_->integerVariable();
4602#endif
4603      for (int i=0;i<numberIntegers;i++) {
4604        CbcSimpleIntegerDynamicPseudoCost * obj =
4605          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(objects[i]) ;
4606        assert (obj&&obj->columnNumber()==integerVariable[i]);
4607        if (info->numberUp_[i]>0) {
4608          if (info->downPseudo_[i]>largest)
4609            largest=info->downPseudo_[i];
4610          if (info->downPseudo_[i]<smallest) 
4611            smallest=info->downPseudo_[i];
4612          if (info->upPseudo_[i]>largest) 
4613            largest=info->upPseudo_[i];
4614          if (info->upPseudo_[i]<smallest)
4615            smallest=info->upPseudo_[i];
4616          obj->updateAfterMini(info->numberDown_[i],
4617                               info->numberDownInfeasible_[i],
4618                               info->downPseudo_[i],
4619                               info->numberUp_[i],
4620                               info->numberUpInfeasible_[i],
4621                               info->upPseudo_[i]);
4622        }
4623      }
4624      //printf("range of costs %g to %g\n",smallest,largest);
4625      simplex->setLogLevel(saveLevel);
4626      numberNodes_ = info->nNodes_;
4627      int numberDo = numberNodes_;
4628      if (whichSolution_>=0)
4629        numberDo--;
4630      if (numberDo>0) {
4631        return 0.5;
4632      } else {
4633        // no solution
4634        return COIN_DBL_MAX; // say infeasible
4635      }
4636    } else {
4637      return -1.0;
4638    }
4639  } else {
4640    return -1.0;
4641  }
4642}
4643
4644// This looks at solution and sets bounds to contain solution
4645void 
4646CbcGeneralDepth::feasibleRegion()
4647{
4648  // Other stuff should have done this
4649}
4650// Redoes data when sequence numbers change
4651void 
4652CbcGeneralDepth::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
4653{
4654}
4655
4656//#define CHECK_PATH
4657#ifdef CHECK_PATH
4658extern const double * debuggerSolution_Z;
4659extern int numberColumns_Z;
4660extern int gotGoodNode_Z;
4661#endif
4662
4663// Creates a branching object
4664CbcBranchingObject * 
4665CbcGeneralDepth::createBranch(int way) 
4666{
4667  int numberDo = numberNodes_;
4668  if (whichSolution_>=0)
4669    numberDo--;
4670  assert (numberDo>0);
4671  // create object
4672  CbcGeneralBranchingObject * branch = new CbcGeneralBranchingObject(model_);
4673  // skip solution
4674  branch->numberSubProblems_ = numberDo;
4675  branch->setNumberBranches(numberDo);
4676  CbcSubProblem * sub = new CbcSubProblem[numberDo];
4677  int iProb=0;
4678  branch->subProblems_ = sub;
4679  branch->numberRows_ = model_->solver()->getNumRows();
4680  int iNode;
4681  OsiSolverInterface * solver = model_->solver();
4682  OsiClpSolverInterface * clpSolver
4683    = dynamic_cast<OsiClpSolverInterface *> (solver);
4684  assert (clpSolver);
4685  ClpSimplex * simplex = clpSolver->getModelPtr();
4686  int numberColumns = simplex->numberColumns();
4687  double * lowerBefore=CoinCopyOfArray(simplex->getColLower(),
4688                                       numberColumns);
4689  double * upperBefore=CoinCopyOfArray(simplex->getColUpper(),
4690                                       numberColumns);
4691  ClpNodeStuff * info = nodeInfo_;
4692  double * weight = new double[numberNodes_];
4693  int * whichNode = new int [numberNodes_];
4694  // Sort
4695  for (iNode=0;iNode<numberNodes_;iNode++) {
4696    if (iNode!=whichSolution_) {
4697      double objectiveValue = info->nodeInfo_[iNode]->objectiveValue();
4698      double sumInfeasibilities = info->nodeInfo_[iNode]->sumInfeasibilities();
4699      int numberInfeasibilities = info->nodeInfo_[iNode]->numberInfeasibilities();
4700      double thisWeight = 0.0;
4701#if 1
4702      // just closest
4703      thisWeight = 1.0e9*numberInfeasibilities;
4704      thisWeight += sumInfeasibilities;
4705      thisWeight += 1.0e-7*objectiveValue;
4706      // Try estimate
4707      thisWeight = info->nodeInfo_[iNode]->estimatedSolution();
4708#else
4709      thisWeight = 1.0e-3*numberInfeasibilities;
4710      thisWeight += 1.0e-5*sumInfeasibilities;
4711      thisWeight += objectiveValue;
4712#endif
4713      whichNode[iProb]=iNode;
4714      weight[iProb++]=thisWeight;
4715    }
4716  }
4717  assert (iProb==numberDo);
4718  CoinSort_2(weight,weight+numberDo,whichNode);
4719  for (iProb=0;iProb<numberDo;iProb++) {
4720    iNode = whichNode[iProb];
4721    ClpNode * node = info->nodeInfo_[iNode];
4722    // move bounds
4723    node->applyNode(simplex,3);
4724    // create subproblem
4725    sub[iProb]=CbcSubProblem(clpSolver,lowerBefore,upperBefore,
4726                             node->statusArray());
4727    sub[iProb].objectiveValue_ = node->objectiveValue();
4728    sub[iProb].sumInfeasibilities_ = node->sumInfeasibilities();
4729    sub[iProb].numberInfeasibilities_ = node->numberInfeasibilities();
4730#ifdef CHECK_PATH
4731    if (simplex->numberColumns()==numberColumns_Z) {
4732      bool onOptimal=true;
4733      const double * columnLower = simplex->columnLower();
4734      const double * columnUpper = simplex->columnUpper();
4735      for (int i=0;i<numberColumns_Z;i++) {
4736        if (iNode==gotGoodNode_Z)
4737          printf("good %d %d %g %g\n",iNode,i,columnLower[i],columnUpper[i]);
4738        if (columnUpper[i]<debuggerSolution_Z[i]||columnLower[i]>debuggerSolution_Z[i]&&simplex->isInteger(i)) {
4739          onOptimal=false;
4740          break;
4741        }
4742      }
4743      if (onOptimal) {
4744        printf("adding to node %x as %d - objs\n",this,iProb);
4745        for (int j=0;j<=iProb;j++)
4746          printf("%d %g\n",j,sub[j].objectiveValue_);
4747      }
4748    }
4749#endif
4750  }
4751  delete [] weight;
4752  delete [] whichNode;
4753  const double * lower = solver->getColLower();
4754  const double * upper = solver->getColUpper();
4755  // restore bounds
4756  for ( int j=0;j<numberColumns;j++) {
4757    if (lowerBefore[j] != lower[j])
4758      solver->setColLower(j,lowerBefore[j]);
4759    if (upperBefore[j] != upper[j])
4760      solver->setColUpper(j,upperBefore[j]);
4761  }
4762  delete [] upperBefore;
4763  delete [] lowerBefore;
4764  return branch;
4765}
4766
4767// Default Constructor
4768CbcGeneralBranchingObject::CbcGeneralBranchingObject()
4769  :CbcBranchingObject(),
4770   subProblems_(NULL),
4771   node_(NULL),
4772   numberSubProblems_(0),
4773   numberRows_(0)
4774{
4775}
4776
4777// Useful constructor
4778CbcGeneralBranchingObject::CbcGeneralBranchingObject (CbcModel * model)
4779  :CbcBranchingObject(model,-1,-1,0.5),
4780   subProblems_(NULL),
4781   node_(NULL),
4782   numberSubProblems_(0),
4783   numberRows_(0)
4784{
4785}
4786
4787// Copy constructor
4788CbcGeneralBranchingObject::CbcGeneralBranchingObject ( const CbcGeneralBranchingObject & rhs) 
4789  :CbcBranchingObject(rhs),
4790   subProblems_(NULL),
4791   node_(rhs.node_),
4792   numberSubProblems_(rhs.numberSubProblems_),
4793   numberRows_(rhs.numberRows_)
4794{
4795  if (numberSubProblems_) {
4796    subProblems_ = new CbcSubProblem[numberSubProblems_];
4797    for (int i=0;i<numberSubProblems_;i++)
4798      subProblems_[i]=rhs.subProblems_[i];
4799  }
4800}
4801
4802// Assignment operator
4803CbcGeneralBranchingObject & 
4804CbcGeneralBranchingObject::operator=( const CbcGeneralBranchingObject& rhs)
4805{
4806  if (this != &rhs) {
4807    CbcBranchingObject::operator=(rhs);
4808    delete [] subProblems_;
4809    numberSubProblems_ = rhs.numberSubProblems_;
4810    numberRows_ = rhs.numberRows_;
4811    if (numberSubProblems_) {
4812      subProblems_ = new CbcSubProblem[numberSubProblems_];
4813      for (int i=0;i<numberSubProblems_;i++)
4814        subProblems_[i]=rhs.subProblems_[i];
4815    } else {
4816      subProblems_ = NULL;
4817    }
4818    node_ = rhs.node_;
4819  }
4820  return *this;
4821}
4822CbcBranchingObject * 
4823CbcGeneralBranchingObject::clone() const
4824{ 
4825  return (new CbcGeneralBranchingObject(*this));
4826}
4827
4828
4829// Destructor
4830CbcGeneralBranchingObject::~CbcGeneralBranchingObject ()
4831{
4832  delete [] subProblems_;
4833}
4834bool doingDoneBranch=false;
4835double
4836CbcGeneralBranchingObject::branch()
4837{
4838  assert (node_);
4839  double cutoff=model_->getCutoff();
4840  bool applied=false;
4841  while (numberBranchesLeft()) {
4842    int which = branchIndex();
4843    decrementNumberBranchesLeft();
4844    CbcSubProblem * thisProb = subProblems_+which;
4845    if (thisProb->objectiveValue_<cutoff) {
4846      //printf("branch %x (sub %x) which now %d\n",this,
4847      //     subProblems_,which);
4848      OsiSolverInterface * solver = model_->solver();
4849      thisProb->apply(solver);
4850      OsiClpSolverInterface * clpSolver
4851        = dynamic_cast<OsiClpSolverInterface *> (solver);
4852      assert (clpSolver);
4853      // Move status to basis
4854      clpSolver->setWarmStart(NULL);
4855      //ClpSimplex * simplex = clpSolver->getModelPtr();
4856      node_->setObjectiveValue(thisProb->objectiveValue_);
4857      node_->setSumInfeasibilities(thisProb->sumInfeasibilities_);
4858      node_->setNumberUnsatisfied(thisProb->numberInfeasibilities_);
4859      applied=true;
4860      doingDoneBranch=true;
4861      break;
4862    } else if (numberBranchesLeft()) {
4863      node_->nodeInfo()->branchedOn() ;
4864    }
4865  }
4866  if (!applied) {
4867    // no good one
4868    node_->setObjectiveValue(cutoff+1.0e20);
4869    node_->setSumInfeasibilities(1.0);
4870    node_->setNumberUnsatisfied(1);
4871  }
4872  return 0.0;
4873}
4874/* Double checks in case node can change its mind!
4875   Can change objective etc */
4876void 
4877CbcGeneralBranchingObject::checkIsCutoff(double cutoff)
4878{
4879  assert (node_);
4880  int first = branchIndex();
4881  int last = first + numberBranchesLeft();
4882  for (int which=first;which<last;which++) {
4883    CbcSubProblem * thisProb = subProblems_+which;
4884    if (thisProb->objectiveValue_<cutoff) {
4885      node_->setObjectiveValue(thisProb->objectiveValue_);
4886      node_->setSumInfeasibilities(thisProb->sumInfeasibilities_);
4887      node_->setNumberUnsatisfied(thisProb->numberInfeasibilities_);
4888      break;
4889    }
4890  }
4891}
4892// Print what would happen 
4893void
4894CbcGeneralBranchingObject::print()
4895{
4896  printf("CbcGeneralObject has %d subproblems\n",numberSubProblems_);
4897}
4898// Fill in current objective etc
4899void 
4900CbcGeneralBranchingObject::state(double & objectiveValue,
4901                                 double & sumInfeasibilities,
4902                                 int & numberUnsatisfied,int which) const
4903{
4904  assert (which>=0&&which<numberSubProblems_);
4905  const CbcSubProblem * thisProb = subProblems_+which;
4906  objectiveValue = thisProb->objectiveValue_;
4907  sumInfeasibilities = thisProb->sumInfeasibilities_;
4908  numberUnsatisfied = thisProb->numberInfeasibilities_;
4909}
4910/** Compare the original object of \c this with the original object of \c
4911    brObj. Assumes that there is an ordering of the original objects.
4912    This method should be invoked only if \c this and brObj are of the same
4913    type.
4914    Return negative/0/positive depending on whether \c this is
4915    smaller/same/larger than the argument.
4916*/
4917int
4918CbcGeneralBranchingObject::compareOriginalObject
4919(const CbcBranchingObject* brObj) const
4920{
4921  throw("must implement");
4922}
4923
4924/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
4925    same type and must have the same original object, but they may have
4926    different feasible regions.
4927    Return the appropriate CbcRangeCompare value (first argument being the
4928    sub/superset if that's the case). In case of overlap (and if \c
4929    replaceIfOverlap is true) replace the current branching object with one
4930    whose feasible region is the overlap.
4931*/
4932CbcRangeCompare
4933CbcGeneralBranchingObject::compareBranchingObject
4934(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
4935{
4936  throw("must implement");
4937}
4938// Default Constructor
4939CbcSubProblem::CbcSubProblem()
4940  : objectiveValue_(0.0),
4941    sumInfeasibilities_(0.0),
4942    variables_(NULL),
4943    newBounds_(NULL),
4944    statusDifference_(NULL),
4945    numberChangedBounds_(0),
4946    numberInfeasibilities_(0)
4947{
4948}
4949
4950// Useful constructor
4951CbcSubProblem::CbcSubProblem (const OsiSolverInterface * solver,
4952                              const double * lastLower,
4953                              const double * lastUpper,
4954                              const unsigned char * status)
4955  : objectiveValue_(0.0),
4956    sumInfeasibilities_(0.0),
4957    variables_(NULL),
4958    newBounds_(NULL),
4959    statusDifference_(NULL),
4960    numberChangedBounds_(0),
4961    numberInfeasibilities_(0)
4962{
4963  const double * lower = solver->getColLower();
4964  const double * upper = solver->getColUpper();
4965
4966  numberChangedBounds_=0;
4967  int numberColumns = solver->getNumCols();
4968  int i;
4969  for (i=0;i<numberColumns;i++) {
4970    if (lower[i]!=lastLower[i]) 
4971      numberChangedBounds_++;
4972    if (upper[i]!=lastUpper[i]) 
4973      numberChangedBounds_++;
4974  }
4975  if (numberChangedBounds_) {
4976    newBounds_ = new double [numberChangedBounds_] ;
4977    variables_ = new int [numberChangedBounds_] ;
4978    numberChangedBounds_=0;
4979    for (i=0;i<numberColumns;i++) {
4980      if (lower[i]!=lastLower[i]) {
4981        variables_[numberChangedBounds_]=i;
4982        newBounds_[numberChangedBounds_++]=lower[i];
4983      }
4984      if (upper[i]!=lastUpper[i]) {
4985        variables_[numberChangedBounds_]=i|0x80000000;
4986        newBounds_[numberChangedBounds_++]=upper[i];
4987      }
4988#ifdef CBC_DEBUG
4989      if (lower[i] != lastLower[i]) {
4990        std::cout
4991          << "lower on " << i << " changed from "
4992          << lastLower[i] << " to " << lower[i] << std::endl ;
4993      }
4994      if (upper[i] != lastUpper[i]) {
4995        std::cout
4996          << "upper on " << i << " changed from "
4997          << lastUpper[i] << " to " << upper[i] << std::endl ;
4998      }
4999#endif
5000    }
5001#ifdef CBC_DEBUG
5002    std::cout << numberChangedBounds_ << " changed bounds." << std::endl ;
5003#endif
5004  }
5005  const OsiClpSolverInterface * clpSolver
5006    = dynamic_cast<const OsiClpSolverInterface *> (solver);
5007  assert (clpSolver);
5008  // Do difference
5009  statusDifference_ = clpSolver->getBasisDiff(status);
5010}
5011
5012// Copy constructor
5013CbcSubProblem::CbcSubProblem ( const CbcSubProblem & rhs) 
5014  : objectiveValue_(rhs.objectiveValue_),
5015    sumInfeasibilities_(rhs.sumInfeasibilities_),
5016    variables_(NULL),
5017    newBounds_(NULL),
5018    statusDifference_(NULL),
5019    numberChangedBounds_(rhs.numberChangedBounds_),
5020    numberInfeasibilities_(rhs.numberInfeasibilities_)
5021{
5022  if (numberChangedBounds_) {
5023    variables_ = CoinCopyOfArray(rhs.variables_,numberChangedBounds_);
5024    newBounds_ = CoinCopyOfArray(rhs.newBounds_,numberChangedBounds_);
5025  }
5026  if (rhs.statusDifference_) {
5027    statusDifference_ = rhs.statusDifference_->clone();
5028  }
5029}
5030
5031// Assignment operator
5032CbcSubProblem & 
5033CbcSubProblem::operator=( const CbcSubProblem& rhs)
5034{
5035  if (this != &rhs) {
5036    delete [] variables_;
5037    delete [] newBounds_;
5038    delete statusDifference_;
5039    objectiveValue_ = rhs.objectiveValue_;
5040    sumInfeasibilities_ = rhs.sumInfeasibilities_;
5041    numberChangedBounds_ = rhs.numberChangedBounds_;
5042    numberInfeasibilities_ = rhs.numberInfeasibilities_;
5043    if (numberChangedBounds_) {
5044      variables_ = CoinCopyOfArray(rhs.variables_,numberChangedBounds_);
5045      newBounds_ = CoinCopyOfArray(rhs.newBounds_,numberChangedBounds_);
5046    } else {
5047      variables_ = NULL;
5048      newBounds_ = NULL;
5049    }
5050    if (rhs.statusDifference_) {
5051      statusDifference_ = rhs.statusDifference_->clone();
5052    } else {
5053      statusDifference_ = NULL;
5054    }
5055  }
5056  return *this;
5057}
5058
5059// Destructor
5060CbcSubProblem::~CbcSubProblem ()
5061{
5062  delete [] variables_;
5063  delete [] newBounds_;
5064  delete statusDifference_;
5065}
5066// Apply subproblem
5067void 
5068CbcSubProblem::apply(OsiSolverInterface * solver)
5069{
5070  int i;
5071  for (i=0;i<numberChangedBounds_;i++) {
5072    int variable = variables_[i];
5073    int k = variable&0x3fffffff;
5074    if ((variable&0x80000000)==0) {
5075      // lower bound changing
5076      //#define CBC_PRINT2
5077#ifdef CBC_PRINT2
5078      if(solver->getColLower()[k]!=newBounds_[i])
5079        printf("lower change for column %d - from %g to %g\n",k,solver->getColLower()[k],newBounds_[i]);
5080#endif
5081#ifndef NDEBUG
5082      if ((variable&0x40000000)==0&&true) {
5083        double oldValue = solver->getColLower()[k];
5084        assert (newBounds_[i]>oldValue-1.0e-8);
5085        if (newBounds_[i]<oldValue+1.0e-8)
5086          printf("bad null lower change for column %d - bound %g\n",k,oldValue);
5087      }
5088#endif
5089      solver->setColLower(k,newBounds_[i]);
5090    } else {
5091      // upper bound changing
5092#ifdef CBC_PRINT2
5093      if(solver->getColUpper()[k]!=newBounds_[i])
5094        printf("upper change for column %d - from %g to %g\n",k,solver->getColUpper()[k],newBounds_[i]);
5095#endif
5096#ifndef NDEBUG
5097      if ((variable&0x40000000)==0&&true) {
5098        double oldValue = solver->getColUpper()[k];
5099        assert (newBounds_[i]<oldValue+1.0e-8);
5100        if (newBounds_[i]>oldValue-1.0e-8)
5101          printf("bad null upper change for column %d - bound %g\n",k,oldValue);
5102      }
5103#endif
5104      solver->setColUpper(k,newBounds_[i]);
5105    }
5106  }
5107  OsiClpSolverInterface * clpSolver
5108    = dynamic_cast<OsiClpSolverInterface *> (solver);
5109  assert (clpSolver);
5110  // Current basis
5111  CoinWarmStartBasis * basis=clpSolver->getPointerToWarmStart();
5112  basis->applyDiff(statusDifference_);
5113  clpSolver->setBasis(*basis);
5114}
5115#endif
5116
5117/** Compare the original object of \c this with the original object of \c
5118    brObj. Assumes that there is an ordering of the original objects.
5119    This method should be invoked only if \c this and brObj are of the same
5120    type.
5121    Return negative/0/positive depending on whether \c this is
5122    smaller/same/larger than the argument.
5123*/
5124int
5125CbcDummyBranchingObject::compareOriginalObject
5126(const CbcBranchingObject* brObj) const
5127{
5128  throw("must implement");
5129}
5130
5131/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
5132    same type and must have the same original object, but they may have
5133    different feasible regions.
5134    Return the appropriate CbcRangeCompare value (first argument being the
5135    sub/superset if that's the case). In case of overlap (and if \c
5136    replaceIfOverlap is true) replace the current branching object with one
5137    whose feasible region is the overlap.
5138*/
5139CbcRangeCompare
5140CbcDummyBranchingObject::compareBranchingObject
5141(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
5142{
5143  throw("must implement");
5144}
Note: See TracBrowser for help on using the repository browser.