source: branches/devel/Cbc/src/CbcBranchActual.cpp @ 477

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

for nonlinear and start moving to OsiTree?
afor n

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 86.7 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 <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "OsiSolverInterface.hpp"
13#include "OsiSolverBranch.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcBranchActual.hpp"
17#include "CoinSort.hpp"
18#include "CoinError.hpp"
19
20// Default Constructor
21CbcClique::CbcClique ()
22  : CbcObject(),
23    numberMembers_(0),
24    numberNonSOSMembers_(0),
25    members_(NULL),
26    type_(NULL),
27    cliqueType_(-1),
28    slack_(-1)
29{
30}
31
32// Useful constructor (which are integer indices)
33CbcClique::CbcClique (CbcModel * model, int cliqueType, int numberMembers,
34           const int * which, const char * type, int identifier,int slack)
35  : CbcObject(model)
36{
37  id_=identifier;
38  numberMembers_=numberMembers;
39  if (numberMembers_) {
40    members_ = new int[numberMembers_];
41    memcpy(members_,which,numberMembers_*sizeof(int));
42    type_ = new char[numberMembers_];
43    if (type) {
44      memcpy(type_,type,numberMembers_*sizeof(char));
45    } else {
46      for (int i=0;i<numberMembers_;i++)
47        type_[i]=1;
48    }
49  } else {
50    members_ = NULL;
51    type_ = NULL;
52  }
53  // Find out how many non sos
54  int i;
55  numberNonSOSMembers_=0;
56  for (i=0;i<numberMembers_;i++)
57    if (!type_[i])
58      numberNonSOSMembers_++;
59  cliqueType_ = cliqueType;
60  slack_ = slack;
61}
62
63// Copy constructor
64CbcClique::CbcClique ( const CbcClique & rhs)
65  :CbcObject(rhs)
66{
67  numberMembers_ = rhs.numberMembers_;
68  numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
69  if (numberMembers_) {
70    members_ = new int[numberMembers_];
71    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
72    type_ = new char[numberMembers_];
73    memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
74  } else {
75    members_ = NULL;
76    type_ = NULL;
77  }
78  cliqueType_ = rhs.cliqueType_;
79  slack_ = rhs.slack_;
80}
81
82// Clone
83CbcObject *
84CbcClique::clone() const
85{
86  return new CbcClique(*this);
87}
88
89// Assignment operator
90CbcClique & 
91CbcClique::operator=( const CbcClique& rhs)
92{
93  if (this!=&rhs) {
94    CbcObject::operator=(rhs);
95    delete [] members_;
96    delete [] type_;
97    numberMembers_ = rhs.numberMembers_;
98    numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
99    if (numberMembers_) {
100      members_ = new int[numberMembers_];
101      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
102      type_ = new char[numberMembers_];
103      memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
104    } else {
105      members_ = NULL;
106      type_ = NULL;
107    }
108    cliqueType_ = rhs.cliqueType_;
109    slack_ = rhs.slack_;
110  }
111  return *this;
112}
113
114// Destructor
115CbcClique::~CbcClique ()
116{
117  delete [] members_;
118  delete [] type_;
119}
120
121// Infeasibility - large is 0.5
122double 
123CbcClique::infeasibility(int & preferredWay) const
124{
125  int numberUnsatis=0, numberFree=0;
126  int j;
127  const int * integer = model_->integerVariable();
128  OsiSolverInterface * solver = model_->solver();
129  const double * solution = model_->testSolution();
130  const double * lower = solver->getColLower();
131  const double * upper = solver->getColUpper();
132  double largestValue=0.0;
133  double integerTolerance = 
134    model_->getDblParam(CbcModel::CbcIntegerTolerance);
135  double * sort = new double[numberMembers_];
136
137  double slackValue=0.0;
138  for (j=0;j<numberMembers_;j++) {
139    int sequence = members_[j];
140    int iColumn = integer[sequence];
141    double value = solution[iColumn];
142    value = CoinMax(value, lower[iColumn]);
143    value = CoinMin(value, upper[iColumn]);
144    double nearest = floor(value+0.5);
145    double distance = fabs(value-nearest);
146    if (distance>integerTolerance) {
147      if (!type_[j])
148        value = 1.0-value; // non SOS
149      // if slack then choose that
150      if (j==slack_&&value>0.05)
151        slackValue = value;
152      largestValue = CoinMax(value,largestValue);
153      sort[numberUnsatis++]=-value;
154    } else if (upper[iColumn]>lower[iColumn]) {
155      numberFree++;
156    }
157  }
158  preferredWay=1;
159  double otherWay = 0.0;
160  if (numberUnsatis) {
161    // sort
162    std::sort(sort,sort+numberUnsatis);
163    for (j=0;j<numberUnsatis;j++) {
164      if ((j&1)!=0)
165        otherWay += -sort[j];
166    }
167    // Need to think more
168    double value = 0.2*numberUnsatis+0.01*(numberMembers_-numberFree);
169    if (fabs(largestValue-0.5)<0.1) {
170      // close to half
171      value +=0.1;
172    }
173    if (slackValue) {
174      // branching on slack
175      value += slackValue;
176    }
177    // scale other way
178    otherWay *= value/(1.0-otherWay);
179    delete [] sort;
180    return value;
181  } else {
182    delete [] sort;
183    return 0.0; // satisfied
184  }
185}
186
187// This looks at solution and sets bounds to contain solution
188void 
189CbcClique::feasibleRegion()
190{
191  int j;
192  const int * integer = model_->integerVariable();
193  OsiSolverInterface * solver = model_->solver();
194  const double * solution = model_->testSolution();
195  const double * lower = solver->getColLower();
196  const double * upper = solver->getColUpper();
197#ifndef NDEBUG
198  double integerTolerance = 
199    model_->getDblParam(CbcModel::CbcIntegerTolerance);
200#endif 
201  for (j=0;j<numberMembers_;j++) {
202    int sequence = members_[j];
203    int iColumn = integer[sequence];
204    double value = solution[iColumn];
205    value = CoinMax(value, lower[iColumn]);
206    value = CoinMin(value, upper[iColumn]);
207    double nearest = floor(value+0.5);
208#ifndef NDEBUG
209    double distance = fabs(value-nearest);
210    assert(distance<=integerTolerance);
211#endif
212    solver->setColLower(iColumn,nearest);
213    solver->setColUpper(iColumn,nearest);
214  }
215}
216// Redoes data when sequence numbers change
217void 
218CbcClique::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
219{
220  model_=model;
221  int n2=0;
222  for (int j=0;j<numberMembers_;j++) {
223    int iColumn = members_[j];
224    int i;
225    for (i=0;i<numberColumns;i++) {
226      if (originalColumns[i]==iColumn)
227        break;
228    }
229    if (i<numberColumns) {
230      members_[n2]=i;
231      type_[n2++]=type_[j];
232    }
233  }
234  if (n2<numberMembers_) {
235    printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
236    numberMembers_=n2;
237  }
238  // Find out how many non sos
239  int i;
240  numberNonSOSMembers_=0;
241  for (i=0;i<numberMembers_;i++)
242    if (!type_[i])
243      numberNonSOSMembers_++;
244}
245
246
247// Creates a branching object
248CbcBranchingObject * 
249CbcClique::createBranch(int way) 
250{
251  int numberUnsatis=0;
252  int j;
253  int nUp=0;
254  int nDown=0;
255  int numberFree=numberMembers_;
256  const int * integer = model_->integerVariable();
257  OsiSolverInterface * solver = model_->solver();
258  const double * solution = model_->testSolution();
259  const double * lower = solver->getColLower();
260  const double * upper = solver->getColUpper();
261  int * upList = new int[numberMembers_];
262  int * downList = new int[numberMembers_];
263  double * sort = new double[numberMembers_];
264  double integerTolerance = 
265      model_->getDblParam(CbcModel::CbcIntegerTolerance);
266
267  double slackValue=0.0;
268  for (j=0;j<numberMembers_;j++) {
269    int sequence = members_[j];
270    int iColumn = integer[sequence];
271    double value = solution[iColumn];
272    value = CoinMax(value, lower[iColumn]);
273    value = CoinMin(value, upper[iColumn]);
274    double nearest = floor(value+0.5);
275    double distance = fabs(value-nearest);
276    if (distance>integerTolerance) {
277      if (!type_[j])
278        value = 1.0-value; // non SOS
279      // if slack then choose that
280      if (j==slack_&&value>0.05)
281        slackValue = value;
282      value = -value; // for sort
283      upList[numberUnsatis]=j;
284      sort[numberUnsatis++]=value;
285    } else if (upper[iColumn]>lower[iColumn]) {
286      upList[--numberFree]=j;
287    }
288  }
289  assert (numberUnsatis);
290  if (!slackValue) {
291    // sort
292    CoinSort_2(sort,sort+numberUnsatis,upList);
293    // put first in up etc
294    int kWay=1;
295    for (j=0;j<numberUnsatis;j++) {
296      if (kWay>0) 
297        upList[nUp++]=upList[j];
298      else
299        downList[nDown++]=upList[j];
300      kWay = -kWay;
301    }
302    for (j=numberFree;j<numberMembers_;j++) {
303      if (kWay>0)
304        upList[nUp++]=upList[j];
305      else
306        downList[nDown++]=upList[j];
307      kWay = -kWay;
308    }
309  } else {
310    // put slack to 0 in first way
311    nUp = 1;
312    upList[0]=slack_;
313    for (j=0;j<numberUnsatis;j++) {
314      downList[nDown++]=upList[j];
315    }
316    for (j=numberFree;j<numberMembers_;j++) {
317      downList[nDown++]=upList[j];
318    }
319  }
320  // create object
321  CbcBranchingObject * branch;
322  if (numberMembers_ <=64)
323     branch = new CbcCliqueBranchingObject(model_,this,way,
324                                         nDown,downList,nUp,upList);
325  else
326    branch = new CbcLongCliqueBranchingObject(model_,this,way,
327                                            nDown,downList,nUp,upList);
328  delete [] upList;
329  delete [] downList;
330  delete [] sort;
331  return branch;
332}
333
334// Default Constructor
335CbcSOS::CbcSOS ()
336  : CbcObject(),
337    members_(NULL),
338    weights_(NULL),
339    numberMembers_(0),
340    sosType_(-1),
341    integerValued_(false)
342{
343}
344
345// Useful constructor (which are indices)
346CbcSOS::CbcSOS (CbcModel * model,  int numberMembers,
347           const int * which, const double * weights, int identifier,int type)
348  : CbcObject(model),
349    numberMembers_(numberMembers),
350    sosType_(type)
351{
352  id_=identifier;
353  integerValued_ = type==1;
354  if (numberMembers_) {
355    members_ = new int[numberMembers_];
356    weights_ = new double[numberMembers_];
357    memcpy(members_,which,numberMembers_*sizeof(int));
358    if (weights) {
359      memcpy(weights_,weights,numberMembers_*sizeof(double));
360    } else {
361      for (int i=0;i<numberMembers_;i++)
362        weights_[i]=i;
363    }
364    // sort so weights increasing
365    CoinSort_2(weights_,weights_+numberMembers_,members_);
366    double last = -COIN_DBL_MAX;
367    int i;
368    for (i=0;i<numberMembers_;i++) {
369      double possible = CoinMax(last+1.0e-10,weights_[i]);
370      weights_[i] = possible;
371      last=possible;
372    }
373  } else {
374    members_ = NULL;
375    weights_ = NULL;
376  }
377  assert (sosType_>0&&sosType_<3);
378}
379
380// Copy constructor
381CbcSOS::CbcSOS ( const CbcSOS & rhs)
382  :CbcObject(rhs)
383{
384  numberMembers_ = rhs.numberMembers_;
385  sosType_ = rhs.sosType_;
386  integerValued_ = rhs.integerValued_;
387  if (numberMembers_) {
388    members_ = new int[numberMembers_];
389    weights_ = new double[numberMembers_];
390    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
391    memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
392  } else {
393    members_ = NULL;
394    weights_ = NULL;
395  }
396}
397
398// Clone
399CbcObject *
400CbcSOS::clone() const
401{
402  return new CbcSOS(*this);
403}
404
405// Assignment operator
406CbcSOS & 
407CbcSOS::operator=( const CbcSOS& rhs)
408{
409  if (this!=&rhs) {
410    CbcObject::operator=(rhs);
411    delete [] members_;
412    delete [] weights_;
413    numberMembers_ = rhs.numberMembers_;
414    sosType_ = rhs.sosType_;
415    integerValued_ = rhs.integerValued_;
416    if (numberMembers_) {
417      members_ = new int[numberMembers_];
418      weights_ = new double[numberMembers_];
419      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
420      memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
421    } else {
422      members_ = NULL;
423      weights_ = NULL;
424    }
425  }
426  return *this;
427}
428
429// Destructor
430CbcSOS::~CbcSOS ()
431{
432  delete [] members_;
433  delete [] weights_;
434}
435
436// Infeasibility - large is 0.5
437double 
438CbcSOS::infeasibility(int & preferredWay) const
439{
440  int j;
441  int firstNonZero=-1;
442  int lastNonZero = -1;
443  OsiSolverInterface * solver = model_->solver();
444  const double * solution = model_->testSolution();
445  //const double * lower = solver->getColLower();
446  const double * upper = solver->getColUpper();
447  //double largestValue=0.0;
448  double integerTolerance = 
449    model_->getDblParam(CbcModel::CbcIntegerTolerance);
450  double weight = 0.0;
451  double sum =0.0;
452
453  // check bounds etc
454  double lastWeight=-1.0e100;
455  for (j=0;j<numberMembers_;j++) {
456    int iColumn = members_[j];
457    if (lastWeight>=weights_[j]-1.0e-7)
458      throw CoinError("Weights too close together in SOS","infeasibility","CbcSOS");
459    double value = CoinMax(0.0,solution[iColumn]);
460    sum += value;
461    if (value>integerTolerance&&upper[iColumn]) {
462      // Possibly due to scaling a fixed variable might slip through
463      if (value>upper[iColumn]) {
464        value=upper[iColumn];
465        // Could change to #ifdef CBC_DEBUG
466#ifndef NDEBUG
467        if (model_->messageHandler()->logLevel()>2)
468          printf("** Variable %d (%d) has value %g and upper bound of %g\n",
469                 iColumn,j,value,upper[iColumn]);
470#endif
471      } 
472      weight += weights_[j]*value;
473      if (firstNonZero<0)
474        firstNonZero=j;
475      lastNonZero=j;
476    }
477  }
478  preferredWay=1;
479  if (lastNonZero-firstNonZero>=sosType_) {
480    // find where to branch
481    assert (sum>0.0);
482    weight /= sum;
483    //int iWhere;
484    //for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++)
485    //if (weight<weights_[iWhere+1])
486    //break;
487    // probably best to use pseudo duals
488    double value = lastNonZero-firstNonZero+1;
489    value *= 0.5/((double) numberMembers_);
490    return value;
491  } else {
492    return 0.0; // satisfied
493  }
494}
495
496// This looks at solution and sets bounds to contain solution
497void 
498CbcSOS::feasibleRegion()
499{
500  int j;
501  int firstNonZero=-1;
502  int lastNonZero = -1;
503  OsiSolverInterface * solver = model_->solver();
504  const double * solution = model_->testSolution();
505  //const double * lower = solver->getColLower();
506  const double * upper = solver->getColUpper();
507  double integerTolerance = 
508    model_->getDblParam(CbcModel::CbcIntegerTolerance);
509  double weight = 0.0;
510  double sum =0.0;
511
512  for (j=0;j<numberMembers_;j++) {
513    int iColumn = members_[j];
514    double value = CoinMax(0.0,solution[iColumn]);
515    sum += value;
516    if (value>integerTolerance&&upper[iColumn]) {
517      weight += weights_[j]*value;
518      if (firstNonZero<0)
519        firstNonZero=j;
520      lastNonZero=j;
521    }
522  }
523  assert (lastNonZero-firstNonZero<sosType_) ;
524  for (j=0;j<firstNonZero;j++) {
525    int iColumn = members_[j];
526    solver->setColUpper(iColumn,0.0);
527  }
528  for (j=lastNonZero+1;j<numberMembers_;j++) {
529    int iColumn = members_[j];
530    solver->setColUpper(iColumn,0.0);
531  }
532}
533// Redoes data when sequence numbers change
534void 
535CbcSOS::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
536{
537  model_=model;
538  int n2=0;
539  for (int j=0;j<numberMembers_;j++) {
540    int iColumn = members_[j];
541    int i;
542    for (i=0;i<numberColumns;i++) {
543      if (originalColumns[i]==iColumn)
544        break;
545    }
546    if (i<numberColumns) {
547      members_[n2]=i;
548      weights_[n2++]=weights_[j];
549    }
550  }
551  if (n2<numberMembers_) {
552    printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
553    numberMembers_=n2;
554  }
555}
556
557
558// Creates a branching object
559CbcBranchingObject * 
560CbcSOS::createBranch(int way) 
561{
562  int j;
563  const double * solution = model_->testSolution();
564  double integerTolerance = 
565      model_->getDblParam(CbcModel::CbcIntegerTolerance);
566  OsiSolverInterface * solver = model_->solver();
567  const double * upper = solver->getColUpper();
568  int firstNonFixed=-1;
569  int lastNonFixed=-1;
570  int firstNonZero=-1;
571  int lastNonZero = -1;
572  double weight = 0.0;
573  double sum =0.0;
574  for (j=0;j<numberMembers_;j++) {
575    int iColumn = members_[j];
576    if (upper[iColumn]) {
577      double value = CoinMax(0.0,solution[iColumn]);
578      sum += value;
579      if (firstNonFixed<0)
580        firstNonFixed=j;
581      lastNonFixed=j;
582      if (value>integerTolerance) {
583        weight += weights_[j]*value;
584        if (firstNonZero<0)
585          firstNonZero=j;
586        lastNonZero=j;
587      }
588    }
589  }
590  assert (lastNonZero-firstNonZero>=sosType_) ;
591  // find where to branch
592  assert (sum>0.0);
593  weight /= sum;
594  int iWhere;
595  double separator=0.0;
596  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
597    if (weight<weights_[iWhere+1])
598      break;
599  if (sosType_==1) {
600    // SOS 1
601    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
602  } else {
603    // SOS 2
604    if (iWhere==firstNonFixed)
605      iWhere++;;
606    if (iWhere==lastNonFixed-1)
607      iWhere = lastNonFixed-2;
608    separator = weights_[iWhere+1];
609  }
610  // create object
611  CbcBranchingObject * branch;
612  branch = new CbcSOSBranchingObject(model_,this,way,separator);
613  return branch;
614}
615
616/* Create an OsiSolverBranch object
617   
618This returns NULL if branch not represented by bound changes
619*/
620OsiSolverBranch * 
621CbcSOS::solverBranch() const
622{
623  int j;
624  const double * solution = model_->testSolution();
625  double integerTolerance = 
626      model_->getDblParam(CbcModel::CbcIntegerTolerance);
627  OsiSolverInterface * solver = model_->solver();
628  const double * upper = solver->getColUpper();
629  int firstNonFixed=-1;
630  int lastNonFixed=-1;
631  int firstNonZero=-1;
632  int lastNonZero = -1;
633  double weight = 0.0;
634  double sum =0.0;
635  double * fix = new double[numberMembers_];
636  int * which = new int[numberMembers_];
637  for (j=0;j<numberMembers_;j++) {
638    int iColumn = members_[j];
639    // fix all on one side or other (even if fixed)
640    fix[j]=0.0;
641    which[j]=iColumn;
642    if (upper[iColumn]) {
643      double value = CoinMax(0.0,solution[iColumn]);
644      sum += value;
645      if (firstNonFixed<0)
646        firstNonFixed=j;
647      lastNonFixed=j;
648      if (value>integerTolerance) {
649        weight += weights_[j]*value;
650        if (firstNonZero<0)
651          firstNonZero=j;
652        lastNonZero=j;
653      }
654    }
655  }
656  assert (lastNonZero-firstNonZero>=sosType_) ;
657  // find where to branch
658  assert (sum>0.0);
659  weight /= sum;
660  // down branch fixes ones above weight to 0
661  int iWhere;
662  int iDownStart=0;
663  int iUpEnd=0;
664  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
665    if (weight<weights_[iWhere+1])
666      break;
667  if (sosType_==1) {
668    // SOS 1
669    iUpEnd=iWhere+1;
670    iDownStart=iUpEnd;
671  } else {
672    // SOS 2
673    if (iWhere==firstNonFixed)
674      iWhere++;;
675    if (iWhere==lastNonFixed-1)
676      iWhere = lastNonFixed-2;
677    iUpEnd=iWhere+1;
678    iDownStart=iUpEnd+1;
679  }
680  //
681  OsiSolverBranch * branch = new OsiSolverBranch();
682  branch->addBranch(-1,0,NULL,NULL,numberMembers_-iDownStart,which+iDownStart,fix);
683  branch->addBranch(1,0,NULL,NULL,iUpEnd,which,fix);
684  delete [] fix;
685  delete [] which;
686  return branch;
687}
688// Construct an OsiSOS object
689OsiSOS * 
690CbcSOS::osiObject(const OsiSolverInterface * solver) const
691{
692  OsiSOS * obj = new OsiSOS(solver,numberMembers_,members_,weights_,sosType_);
693  obj->setPriority(priority());
694  return obj;
695}
696
697/** Default Constructor
698
699  Equivalent to an unspecified binary variable.
700*/
701CbcSimpleInteger::CbcSimpleInteger ()
702  : CbcObject(),
703    originalLower_(0.0),
704    originalUpper_(1.0),
705    breakEven_(0.5),
706    columnNumber_(-1),
707    preferredWay_(0)
708{
709}
710
711/** Useful constructor
712
713  Loads actual upper & lower bounds for the specified variable.
714*/
715CbcSimpleInteger::CbcSimpleInteger ( CbcModel * model, int iColumn, double breakEven)
716  : CbcObject(model)
717{
718  columnNumber_ = iColumn ;
719  originalLower_ = model->solver()->getColLower()[columnNumber_] ;
720  originalUpper_ = model->solver()->getColUpper()[columnNumber_] ;
721  breakEven_ = breakEven;
722  assert (breakEven_>0.0&&breakEven_<1.0);
723  preferredWay_ = 0;
724}
725
726
727// Copy constructor
728CbcSimpleInteger::CbcSimpleInteger ( const CbcSimpleInteger & rhs)
729  :CbcObject(rhs)
730
731{
732  columnNumber_ = rhs.columnNumber_;
733  originalLower_ = rhs.originalLower_;
734  originalUpper_ = rhs.originalUpper_;
735  breakEven_ = rhs.breakEven_;
736  preferredWay_ = rhs.preferredWay_;
737}
738
739// Clone
740CbcObject *
741CbcSimpleInteger::clone() const
742{
743  return new CbcSimpleInteger(*this);
744}
745
746// Assignment operator
747CbcSimpleInteger & 
748CbcSimpleInteger::operator=( const CbcSimpleInteger& rhs)
749{
750  if (this!=&rhs) {
751    CbcObject::operator=(rhs);
752    columnNumber_ = rhs.columnNumber_;
753    originalLower_ = rhs.originalLower_;
754    originalUpper_ = rhs.originalUpper_;
755    breakEven_ = rhs.breakEven_;
756    preferredWay_ = rhs.preferredWay_;
757  }
758  return *this;
759}
760
761// Destructor
762CbcSimpleInteger::~CbcSimpleInteger ()
763{
764}
765// Construct an OsiSimpleInteger object
766OsiSimpleInteger * 
767CbcSimpleInteger::osiObject() const
768{
769  OsiSimpleInteger * obj = new OsiSimpleInteger(columnNumber_,
770                                                originalLower_,originalUpper_);
771  obj->setPriority(priority());
772  return obj;
773}
774
775double
776CbcSimpleInteger::infeasibility(const OsiSolverInterface * solver, const OsiBranchingInformation * info,
777                         int & preferredWay) const
778{
779  double value = info->solution_[columnNumber_];
780  value = CoinMax(value, info->lower_[columnNumber_]);
781  value = CoinMin(value, info->upper_[columnNumber_]);
782  double nearest = floor(value+(1.0-breakEven_));
783  assert (breakEven_>0.0&&breakEven_<1.0);
784  if (nearest>value) 
785    preferredWay=1;
786  else
787    preferredWay=-1;
788  if (preferredWay_)
789    preferredWay=preferredWay_;
790  double weight = fabs(value-nearest);
791  // normalize so weight is 0.5 at break even
792  if (nearest<value)
793    weight = (0.5/breakEven_)*weight;
794  else
795    weight = (0.5/(1.0-breakEven_))*weight;
796  if (fabs(value-nearest)<=info->integerTolerance_) 
797    return 0.0;
798  else
799    return weight;
800}
801double 
802CbcSimpleInteger::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const 
803{
804  double value = info->solution_[columnNumber_];
805  double newValue = CoinMax(value, info->lower_[columnNumber_]);
806  newValue = CoinMin(newValue, info->upper_[columnNumber_]);
807  newValue = floor(newValue+0.5);
808  solver->setColLower(columnNumber_,newValue);
809  solver->setColUpper(columnNumber_,newValue);
810  return fabs(value-newValue);
811}
812
813/* Create an OsiSolverBranch object
814   
815This returns NULL if branch not represented by bound changes
816*/
817OsiSolverBranch * 
818CbcSimpleInteger::solverBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
819{
820  double value = info->solution_[columnNumber_];
821  value = CoinMax(value, info->lower_[columnNumber_]);
822  value = CoinMin(value, info->upper_[columnNumber_]);
823  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
824#ifndef NDEBUG
825  double nearest = floor(value+0.5);
826  assert (fabs(value-nearest)>info->integerTolerance_);
827#endif
828  OsiSolverBranch * branch = new OsiSolverBranch();
829  branch->addBranch(columnNumber_,value);
830  return branch;
831}
832// Creates a branching object
833CbcBranchingObject * 
834CbcSimpleInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) 
835{
836  double value = info->solution_[columnNumber_];
837  value = CoinMax(value, info->lower_[columnNumber_]);
838  value = CoinMin(value, info->upper_[columnNumber_]);
839  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
840  if (!info->hotstartSolution_) {
841#ifndef NDEBUG
842    double nearest = floor(value+0.5);
843    assert (fabs(value-nearest)>info->integerTolerance_);
844#endif
845  } else {
846    double targetValue = info->hotstartSolution_[columnNumber_];
847    if (way>0)
848      value = targetValue-0.1;
849    else
850      value = targetValue+0.1;
851  }
852  CbcBranchingObject * branch = new CbcIntegerBranchingObject(model_,columnNumber_,way,
853                                             value);
854  branch->setOriginalObject(this);
855  return branch;
856}
857/* Column number if single column object -1 otherwise,
858   so returns >= 0
859   Used by heuristics
860*/
861int 
862CbcSimpleInteger::columnNumber() const
863{
864  return columnNumber_;
865}
866/* Reset variable bounds to their original values.
867 
868    Bounds may be tightened, so it may be good to be able to set this info in object.
869*/
870void 
871CbcSimpleInteger::resetBounds(const OsiSolverInterface * solver) 
872{
873  originalLower_ = solver->getColLower()[columnNumber_] ;
874  originalUpper_ = solver->getColUpper()[columnNumber_] ;
875}
876
877// Infeasibility - large is 0.5
878double 
879CbcSimpleInteger::infeasibility(int & preferredWay) const
880{
881  OsiBranchingInformation info(model_->solver(),model_->normalSolver());
882  return infeasibility(model_->solver(),&info,preferredWay);
883}
884
885// This looks at solution and sets bounds to contain solution
886/** More precisely: it first forces the variable within the existing
887    bounds, and then tightens the bounds to fix the variable at the
888    nearest integer value.
889*/
890void 
891CbcSimpleInteger::feasibleRegion()
892{
893  abort();
894}
895CbcBranchingObject * 
896CbcSimpleInteger::createBranch( int way) 
897{
898  abort();
899  return NULL;
900}
901// Default Constructor
902CbcIntegerBranchingObject::CbcIntegerBranchingObject()
903  :CbcBranchingObject()
904{
905  down_[0] = 0.0;
906  down_[1] = 0.0;
907  up_[0] = 0.0;
908  up_[1] = 0.0;
909}
910
911// Useful constructor
912CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
913                                                      int variable, int way , double value)
914  :CbcBranchingObject(model,variable,way,value)
915{
916  int iColumn = variable;
917  down_[0] = model_->solver()->getColLower()[iColumn];
918  down_[1] = floor(value_);
919  up_[0] = ceil(value_);
920  up_[1] = model->getColUpper()[iColumn];
921}
922// Useful constructor for fixing
923CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
924                                                      int variable, int way,
925                                                      double lowerValue, 
926                                                      double upperValue)
927  :CbcBranchingObject(model,variable,way,lowerValue)
928{
929  setNumberBranchesLeft(1);
930  down_[0] = lowerValue;
931  down_[1] = upperValue;
932  up_[0] = lowerValue;
933  up_[1] = upperValue;
934}
935 
936
937// Copy constructor
938CbcIntegerBranchingObject::CbcIntegerBranchingObject ( const CbcIntegerBranchingObject & rhs) :CbcBranchingObject(rhs)
939{
940  down_[0] = rhs.down_[0];
941  down_[1] = rhs.down_[1];
942  up_[0] = rhs.up_[0];
943  up_[1] = rhs.up_[1];
944}
945
946// Assignment operator
947CbcIntegerBranchingObject & 
948CbcIntegerBranchingObject::operator=( const CbcIntegerBranchingObject& rhs)
949{
950  if (this != &rhs) {
951    CbcBranchingObject::operator=(rhs);
952    down_[0] = rhs.down_[0];
953    down_[1] = rhs.down_[1];
954    up_[0] = rhs.up_[0];
955    up_[1] = rhs.up_[1];
956  }
957  return *this;
958}
959CbcBranchingObject * 
960CbcIntegerBranchingObject::clone() const
961{ 
962  return (new CbcIntegerBranchingObject(*this));
963}
964
965
966// Destructor
967CbcIntegerBranchingObject::~CbcIntegerBranchingObject ()
968{
969}
970
971/*
972  Perform a branch by adjusting the bounds of the specified variable. Note
973  that each arm of the branch advances the object to the next arm by
974  advancing the value of way_.
975
976  Providing new values for the variable's lower and upper bounds for each
977  branching direction gives a little bit of additional flexibility and will
978  be easily extensible to multi-way branching.
979  Returns change in guessed objective on next branch
980*/
981double
982CbcIntegerBranchingObject::branch()
983{
984  decrementNumberBranchesLeft();
985  int iColumn = originalCbcObject_->columnNumber();
986  assert (variable_==iColumn);
987  double olb,oub ;
988  olb = model_->solver()->getColLower()[iColumn] ;
989  oub = model_->solver()->getColUpper()[iColumn] ;
990  if (way_<0) {
991#ifdef CBC_DEBUG
992  { double olb,oub ;
993    olb = model_->solver()->getColLower()[iColumn] ;
994    oub = model_->solver()->getColUpper()[iColumn] ;
995    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
996           iColumn,olb,oub,down_[0],down_[1]) ; }
997#endif
998    model_->solver()->setColLower(iColumn,down_[0]);
999    model_->solver()->setColUpper(iColumn,down_[1]);
1000    way_=1;
1001  } else {
1002#ifdef CBC_DEBUG
1003  { double olb,oub ;
1004    olb = model_->solver()->getColLower()[iColumn] ;
1005    oub = model_->solver()->getColUpper()[iColumn] ;
1006    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
1007           iColumn,olb,oub,up_[0],up_[1]) ; }
1008#endif
1009    model_->solver()->setColLower(iColumn,up_[0]);
1010    model_->solver()->setColUpper(iColumn,up_[1]);
1011    way_=-1;      // Swap direction
1012  }
1013  double nlb = model_->solver()->getColLower()[iColumn];
1014  if (nlb<olb) {
1015#ifndef NDEBUG
1016    printf("bad lb change for column %d from %g to %g\n",iColumn,olb,nlb);
1017#endif
1018    model_->solver()->setColLower(iColumn,olb);
1019  }
1020  double nub = model_->solver()->getColUpper()[iColumn];
1021  if (nub>oub) {
1022#ifndef NDEBUG
1023    printf("bad ub change for column %d from %g to %g\n",iColumn,oub,nub);
1024#endif
1025    model_->solver()->setColUpper(iColumn,oub);
1026  }
1027#ifndef NDEBUG
1028  if (nlb<olb+1.0e-8&&nub>oub-1.0e-8)
1029    printf("bad null change for column %d - bounds %g,%g\n",iColumn,olb,oub);
1030#endif
1031  return 0.0;
1032}
1033// Print what would happen 
1034void
1035CbcIntegerBranchingObject::print()
1036{
1037  int iColumn = originalCbcObject_->columnNumber();
1038  assert (variable_==iColumn);
1039  if (way_<0) {
1040  { double olb,oub ;
1041    olb = model_->solver()->getColLower()[iColumn] ;
1042    oub = model_->solver()->getColUpper()[iColumn] ;
1043    printf("CbcInteger would branch down on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1044           iColumn,variable_,olb,oub,down_[0],down_[1]) ; }
1045  } else {
1046  { double olb,oub ;
1047    olb = model_->solver()->getColLower()[iColumn] ;
1048    oub = model_->solver()->getColUpper()[iColumn] ;
1049    printf("CbcInteger would branch up on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1050           iColumn,variable_,olb,oub,up_[0],up_[1]) ; }
1051  }
1052}
1053
1054
1055/** Default Constructor
1056
1057  Equivalent to an unspecified binary variable.
1058*/
1059CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ()
1060  : CbcSimpleInteger(),
1061    downPseudoCost_(1.0e-5),
1062    upPseudoCost_(1.0e-5),
1063    upDownSeparator_(-1.0),
1064    method_(0)
1065{
1066}
1067
1068/** Useful constructor
1069
1070  Loads actual upper & lower bounds for the specified variable.
1071*/
1072CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1073                                    int iColumn, double breakEven)
1074  : CbcSimpleInteger(model,iColumn,breakEven)
1075{
1076  const double * cost = model->getObjCoefficients();
1077  double costValue = CoinMax(1.0e-5,fabs(cost[iColumn]));
1078  // treat as if will cost what it says up
1079  upPseudoCost_=costValue;
1080  // and balance at breakeven
1081  downPseudoCost_=((1.0-breakEven_)*upPseudoCost_)/breakEven_;
1082  upDownSeparator_ = -1.0;
1083  method_=0;
1084}
1085
1086/** Useful constructor
1087
1088  Loads actual upper & lower bounds for the specified variable.
1089*/
1090CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1091                                    int iColumn, double downPseudoCost,
1092                                                        double upPseudoCost)
1093  : CbcSimpleInteger(model,iColumn)
1094{
1095  downPseudoCost_ = CoinMax(1.0e-10,downPseudoCost);
1096  upPseudoCost_ = CoinMax(1.0e-10,upPseudoCost);
1097  breakEven_ = upPseudoCost_/(upPseudoCost_+downPseudoCost_);
1098  upDownSeparator_ = -1.0;
1099  method_=0;
1100}
1101// Useful constructor - passed and model index and pseudo costs
1102CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model, int dummy,int iColumn, 
1103                                                        double downPseudoCost, double upPseudoCost)
1104{
1105  CbcSimpleIntegerPseudoCost(model,iColumn,downPseudoCost,upPseudoCost);
1106}
1107
1108// Copy constructor
1109CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ( const CbcSimpleIntegerPseudoCost & rhs)
1110  :CbcSimpleInteger(rhs),
1111   downPseudoCost_(rhs.downPseudoCost_),
1112   upPseudoCost_(rhs.upPseudoCost_),
1113   upDownSeparator_(rhs.upDownSeparator_),
1114   method_(rhs.method_)
1115
1116{
1117}
1118
1119// Clone
1120CbcObject *
1121CbcSimpleIntegerPseudoCost::clone() const
1122{
1123  return new CbcSimpleIntegerPseudoCost(*this);
1124}
1125
1126// Assignment operator
1127CbcSimpleIntegerPseudoCost & 
1128CbcSimpleIntegerPseudoCost::operator=( const CbcSimpleIntegerPseudoCost& rhs)
1129{
1130  if (this!=&rhs) {
1131    CbcSimpleInteger::operator=(rhs);
1132    downPseudoCost_=rhs.downPseudoCost_;
1133    upPseudoCost_=rhs.upPseudoCost_;
1134    upDownSeparator_=rhs.upDownSeparator_;
1135    method_=rhs.method_;
1136  }
1137  return *this;
1138}
1139
1140// Destructor
1141CbcSimpleIntegerPseudoCost::~CbcSimpleIntegerPseudoCost ()
1142{
1143}
1144// Creates a branching object
1145CbcBranchingObject * 
1146CbcSimpleIntegerPseudoCost::createBranch(int way) 
1147{
1148  OsiSolverInterface * solver = model_->solver();
1149  const double * solution = model_->testSolution();
1150  const double * lower = solver->getColLower();
1151  const double * upper = solver->getColUpper();
1152  double value = solution[columnNumber_];
1153  value = CoinMax(value, lower[columnNumber_]);
1154  value = CoinMin(value, upper[columnNumber_]);
1155#ifndef NDEBUG
1156  double nearest = floor(value+0.5);
1157  double integerTolerance = 
1158    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1159  assert (upper[columnNumber_]>lower[columnNumber_]);
1160#endif
1161  if (!model_->hotstartSolution()) {
1162    assert (fabs(value-nearest)>integerTolerance);
1163  } else {
1164    const double * hotstartSolution = model_->hotstartSolution();
1165    double targetValue = hotstartSolution[columnNumber_];
1166    if (way>0)
1167      value = targetValue-0.1;
1168    else
1169      value = targetValue+0.1;
1170  }
1171  CbcIntegerPseudoCostBranchingObject * newObject = 
1172    new CbcIntegerPseudoCostBranchingObject(model_,columnNumber_,way,
1173                                            value);
1174  double up =  upPseudoCost_*(ceil(value)-value);
1175  double down =  downPseudoCost_*(value-floor(value));
1176  double changeInGuessed=up-down;
1177  if (way>0)
1178    changeInGuessed = - changeInGuessed;
1179  changeInGuessed=CoinMax(0.0,changeInGuessed);
1180  //if (way>0)
1181  //changeInGuessed += 1.0e8; // bias to stay up
1182  newObject->setChangeInGuessed(changeInGuessed);
1183  newObject->setOriginalObject(this);
1184  return newObject;
1185}
1186// Infeasibility - large is 0.5
1187double 
1188CbcSimpleIntegerPseudoCost::infeasibility(int & preferredWay) const
1189{
1190  OsiSolverInterface * solver = model_->solver();
1191  const double * solution = model_->testSolution();
1192  const double * lower = solver->getColLower();
1193  const double * upper = solver->getColUpper();
1194  if (upper[columnNumber_]==lower[columnNumber_]) {
1195    // fixed
1196    preferredWay=1;
1197    return 0.0;
1198  }
1199  double value = solution[columnNumber_];
1200  value = CoinMax(value, lower[columnNumber_]);
1201  value = CoinMin(value, upper[columnNumber_]);
1202  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
1203    solution[columnNumber_],upper[columnNumber_]);*/
1204  double nearest = floor(value+0.5);
1205  double integerTolerance = 
1206    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1207  double below = floor(value+integerTolerance);
1208  double above = below+1.0;
1209  if (above>upper[columnNumber_]) {
1210    above=below;
1211    below = above -1;
1212  }
1213  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1214  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1215  if (downCost>=upCost)
1216    preferredWay=1;
1217  else
1218    preferredWay=-1;
1219  // See if up down choice set
1220  if (upDownSeparator_>0.0) {
1221    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
1222  }
1223  if (preferredWay_)
1224    preferredWay=preferredWay_;
1225  if (fabs(value-nearest)<=integerTolerance) {
1226    return 0.0;
1227  } else {
1228    // can't get at model so 1,2 don't make sense
1229    assert(method_<1||method_>2);
1230    if (!method_)
1231      return CoinMin(downCost,upCost);
1232    else
1233      return CoinMax(downCost,upCost);
1234  }
1235}
1236
1237// Return "up" estimate
1238double 
1239CbcSimpleIntegerPseudoCost::upEstimate() const
1240{
1241  OsiSolverInterface * solver = model_->solver();
1242  const double * solution = model_->testSolution();
1243  const double * lower = solver->getColLower();
1244  const double * upper = solver->getColUpper();
1245  double value = solution[columnNumber_];
1246  value = CoinMax(value, lower[columnNumber_]);
1247  value = CoinMin(value, upper[columnNumber_]);
1248  if (upper[columnNumber_]==lower[columnNumber_]) {
1249    // fixed
1250    return 0.0;
1251  }
1252  double integerTolerance = 
1253    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1254  double below = floor(value+integerTolerance);
1255  double above = below+1.0;
1256  if (above>upper[columnNumber_]) {
1257    above=below;
1258    below = above -1;
1259  }
1260  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1261  return upCost;
1262}
1263// Return "down" estimate
1264double 
1265CbcSimpleIntegerPseudoCost::downEstimate() const
1266{
1267  OsiSolverInterface * solver = model_->solver();
1268  const double * solution = model_->testSolution();
1269  const double * lower = solver->getColLower();
1270  const double * upper = solver->getColUpper();
1271  double value = solution[columnNumber_];
1272  value = CoinMax(value, lower[columnNumber_]);
1273  value = CoinMin(value, upper[columnNumber_]);
1274  if (upper[columnNumber_]==lower[columnNumber_]) {
1275    // fixed
1276    return 0.0;
1277  }
1278  double integerTolerance = 
1279    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1280  double below = floor(value+integerTolerance);
1281  double above = below+1.0;
1282  if (above>upper[columnNumber_]) {
1283    above=below;
1284    below = above -1;
1285  }
1286  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1287  return downCost;
1288}
1289
1290// Default Constructor
1291CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1292  :CbcIntegerBranchingObject()
1293{
1294  changeInGuessed_=1.0e-5;
1295}
1296
1297// Useful constructor
1298CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1299                                                      int variable, int way , double value)
1300  :CbcIntegerBranchingObject(model,variable,way,value)
1301{
1302}
1303// Useful constructor for fixing
1304CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1305                                                      int variable, int way,
1306                                                      double lowerValue, 
1307                                                      double upperValue)
1308  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
1309{
1310  changeInGuessed_=1.0e100;
1311}
1312 
1313
1314// Copy constructor
1315CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject ( 
1316                                 const CbcIntegerPseudoCostBranchingObject & rhs)
1317  :CbcIntegerBranchingObject(rhs)
1318{
1319  changeInGuessed_ = rhs.changeInGuessed_;
1320}
1321
1322// Assignment operator
1323CbcIntegerPseudoCostBranchingObject & 
1324CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject& rhs)
1325{
1326  if (this != &rhs) {
1327    CbcIntegerBranchingObject::operator=(rhs);
1328    changeInGuessed_ = rhs.changeInGuessed_;
1329  }
1330  return *this;
1331}
1332CbcBranchingObject * 
1333CbcIntegerPseudoCostBranchingObject::clone() const
1334{ 
1335  return (new CbcIntegerPseudoCostBranchingObject(*this));
1336}
1337
1338
1339// Destructor
1340CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1341{
1342}
1343
1344/*
1345  Perform a branch by adjusting the bounds of the specified variable. Note
1346  that each arm of the branch advances the object to the next arm by
1347  advancing the value of way_.
1348
1349  Providing new values for the variable's lower and upper bounds for each
1350  branching direction gives a little bit of additional flexibility and will
1351  be easily extensible to multi-way branching.
1352  Returns change in guessed objective on next branch
1353*/
1354double
1355CbcIntegerPseudoCostBranchingObject::branch()
1356{
1357  CbcIntegerBranchingObject::branch();
1358  return changeInGuessed_;
1359}
1360
1361
1362// Default Constructor
1363CbcCliqueBranchingObject::CbcCliqueBranchingObject()
1364  :CbcBranchingObject()
1365{
1366  clique_ = NULL;
1367  downMask_[0]=0;
1368  downMask_[1]=0;
1369  upMask_[0]=0;
1370  upMask_[1]=0;
1371}
1372
1373// Useful constructor
1374CbcCliqueBranchingObject::CbcCliqueBranchingObject (CbcModel * model,
1375                                                    const CbcClique * clique,
1376                                                    int way ,
1377                                                    int numberOnDownSide, const int * down,
1378                                                    int numberOnUpSide, const int * up)
1379  :CbcBranchingObject(model,clique->id(),way,0.5)
1380{
1381  clique_ = clique;
1382  downMask_[0]=0;
1383  downMask_[1]=0;
1384  upMask_[0]=0;
1385  upMask_[1]=0;
1386  int i;
1387  for (i=0;i<numberOnDownSide;i++) {
1388    int sequence = down[i];
1389    int iWord = sequence>>5;
1390    int iBit = sequence - 32*iWord;
1391    unsigned int k = 1<<iBit;
1392    downMask_[iWord] |= k;
1393  }
1394  for (i=0;i<numberOnUpSide;i++) {
1395    int sequence = up[i];
1396    int iWord = sequence>>5;
1397    int iBit = sequence - 32*iWord;
1398    unsigned int k = 1<<iBit;
1399    upMask_[iWord] |= k;
1400  }
1401}
1402
1403// Copy constructor
1404CbcCliqueBranchingObject::CbcCliqueBranchingObject ( const CbcCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1405{
1406  clique_=rhs.clique_;
1407  downMask_[0]=rhs.downMask_[0];
1408  downMask_[1]=rhs.downMask_[1];
1409  upMask_[0]=rhs.upMask_[0];
1410  upMask_[1]=rhs.upMask_[1];
1411}
1412
1413// Assignment operator
1414CbcCliqueBranchingObject & 
1415CbcCliqueBranchingObject::operator=( const CbcCliqueBranchingObject& rhs)
1416{
1417  if (this != &rhs) {
1418    CbcBranchingObject::operator=(rhs);
1419    clique_=rhs.clique_;
1420    downMask_[0]=rhs.downMask_[0];
1421    downMask_[1]=rhs.downMask_[1];
1422    upMask_[0]=rhs.upMask_[0];
1423    upMask_[1]=rhs.upMask_[1];
1424  }
1425  return *this;
1426}
1427CbcBranchingObject * 
1428CbcCliqueBranchingObject::clone() const
1429{ 
1430  return (new CbcCliqueBranchingObject(*this));
1431}
1432
1433
1434// Destructor
1435CbcCliqueBranchingObject::~CbcCliqueBranchingObject ()
1436{
1437}
1438double
1439CbcCliqueBranchingObject::branch()
1440{
1441  decrementNumberBranchesLeft();
1442  int iWord;
1443  int numberMembers = clique_->numberMembers();
1444  const int * which = clique_->members();
1445  const int * integerVariables = model_->integerVariable();
1446  int numberWords=(numberMembers+31)>>5;
1447  // *** for way - up means fix all those in down section
1448  if (way_<0) {
1449#ifdef FULL_PRINT
1450    printf("Down Fix ");
1451#endif
1452    for (iWord=0;iWord<numberWords;iWord++) {
1453      int i;
1454      for (i=0;i<32;i++) {
1455        unsigned int k = 1<<i;
1456        if ((upMask_[iWord]&k)!=0) {
1457          int iColumn = which[i+32*iWord];
1458#ifdef FULL_PRINT
1459          printf("%d ",i+32*iWord);
1460#endif
1461          // fix weak way
1462          if (clique_->type(i+32*iWord))
1463            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1464          else
1465            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1466        }
1467      }
1468    }
1469    way_=1;       // Swap direction
1470  } else {
1471#ifdef FULL_PRINT
1472    printf("Up Fix ");
1473#endif
1474    for (iWord=0;iWord<numberWords;iWord++) {
1475      int i;
1476      for (i=0;i<32;i++) {
1477        unsigned int k = 1<<i;
1478        if ((downMask_[iWord]&k)!=0) {
1479          int iColumn = which[i+32*iWord];
1480#ifdef FULL_PRINT
1481          printf("%d ",i+32*iWord);
1482#endif
1483          // fix weak way
1484          if (clique_->type(i+32*iWord))
1485            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1486          else
1487            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1488        }
1489      }
1490    }
1491    way_=-1;      // Swap direction
1492  }
1493#ifdef FULL_PRINT
1494  printf("\n");
1495#endif
1496  return 0.0;
1497}
1498// Print what would happen 
1499void
1500CbcCliqueBranchingObject::print()
1501{
1502  int iWord;
1503  int numberMembers = clique_->numberMembers();
1504  const int * which = clique_->members();
1505  const int * integerVariables = model_->integerVariable();
1506  int numberWords=(numberMembers+31)>>5;
1507  // *** for way - up means fix all those in down section
1508  if (way_<0) {
1509    printf("Clique - Down Fix ");
1510    for (iWord=0;iWord<numberWords;iWord++) {
1511      int i;
1512      for (i=0;i<32;i++) {
1513        unsigned int k = 1<<i;
1514        if ((upMask_[iWord]&k)!=0) {
1515          int iColumn = which[i+32*iWord];
1516          printf("%d ",integerVariables[iColumn]);
1517        }
1518      }
1519    }
1520  } else {
1521    printf("Clique - Up Fix ");
1522    for (iWord=0;iWord<numberWords;iWord++) {
1523      int i;
1524      for (i=0;i<32;i++) {
1525        unsigned int k = 1<<i;
1526        if ((downMask_[iWord]&k)!=0) {
1527          int iColumn = which[i+32*iWord];
1528          printf("%d ",integerVariables[iColumn]);
1529        }
1530      }
1531    }
1532  }
1533  printf("\n");
1534}
1535 
1536// Default Constructor
1537CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject()
1538  :CbcBranchingObject()
1539{
1540  clique_=NULL;
1541  downMask_=NULL;
1542  upMask_=NULL;
1543}
1544
1545// Useful constructor
1546CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject (CbcModel * model,
1547                                                            const CbcClique * clique, 
1548                                                            int way ,
1549                                                    int numberOnDownSide, const int * down,
1550                                                    int numberOnUpSide, const int * up)
1551  :CbcBranchingObject(model,clique->id(),way,0.5)
1552{
1553  clique_ = clique;
1554  int numberMembers = clique_->numberMembers();
1555  int numberWords=(numberMembers+31)>>5;
1556  downMask_ = new unsigned int [numberWords];
1557  upMask_ = new unsigned int [numberWords];
1558  memset(downMask_,0,numberWords*sizeof(unsigned int));
1559  memset(upMask_,0,numberWords*sizeof(unsigned int));
1560  int i;
1561  for (i=0;i<numberOnDownSide;i++) {
1562    int sequence = down[i];
1563    int iWord = sequence>>5;
1564    int iBit = sequence - 32*iWord;
1565    unsigned int k = 1<<iBit;
1566    downMask_[iWord] |= k;
1567  }
1568  for (i=0;i<numberOnUpSide;i++) {
1569    int sequence = up[i];
1570    int iWord = sequence>>5;
1571    int iBit = sequence - 32*iWord;
1572    unsigned int k = 1<<iBit;
1573    upMask_[iWord] |= k;
1574  }
1575}
1576
1577// Copy constructor
1578CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject ( const CbcLongCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1579{
1580  clique_=rhs.clique_;
1581  if (rhs.downMask_) {
1582    int numberMembers = clique_->numberMembers();
1583    int numberWords=(numberMembers+31)>>5;
1584    downMask_ = new unsigned int [numberWords];
1585    memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
1586    upMask_ = new unsigned int [numberWords];
1587    memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
1588  } else {
1589    downMask_=NULL;
1590    upMask_=NULL;
1591  }   
1592}
1593
1594// Assignment operator
1595CbcLongCliqueBranchingObject & 
1596CbcLongCliqueBranchingObject::operator=( const CbcLongCliqueBranchingObject& rhs)
1597{
1598  if (this != &rhs) {
1599    CbcBranchingObject::operator=(rhs);
1600    clique_=rhs.clique_;
1601    delete [] downMask_;
1602    delete [] upMask_;
1603    if (rhs.downMask_) {
1604      int numberMembers = clique_->numberMembers();
1605      int numberWords=(numberMembers+31)>>5;
1606      downMask_ = new unsigned int [numberWords];
1607      memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
1608      upMask_ = new unsigned int [numberWords];
1609      memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
1610    } else {
1611      downMask_=NULL;
1612      upMask_=NULL;
1613    }   
1614  }
1615  return *this;
1616}
1617CbcBranchingObject * 
1618CbcLongCliqueBranchingObject::clone() const
1619{ 
1620  return (new CbcLongCliqueBranchingObject(*this));
1621}
1622
1623
1624// Destructor
1625CbcLongCliqueBranchingObject::~CbcLongCliqueBranchingObject ()
1626{
1627  delete [] downMask_;
1628  delete [] upMask_;
1629}
1630double
1631CbcLongCliqueBranchingObject::branch()
1632{
1633  decrementNumberBranchesLeft();
1634  int iWord;
1635  int numberMembers = clique_->numberMembers();
1636  const int * which = clique_->members();
1637  const int * integerVariables = model_->integerVariable();
1638  int numberWords=(numberMembers+31)>>5;
1639  // *** for way - up means fix all those in down section
1640  if (way_<0) {
1641#ifdef FULL_PRINT
1642    printf("Down Fix ");
1643#endif
1644    for (iWord=0;iWord<numberWords;iWord++) {
1645      int i;
1646      for (i=0;i<32;i++) {
1647        unsigned int k = 1<<i;
1648        if ((upMask_[iWord]&k)!=0) {
1649          int iColumn = which[i+32*iWord];
1650#ifdef FULL_PRINT
1651          printf("%d ",i+32*iWord);
1652#endif
1653          // fix weak way
1654          if (clique_->type(i+32*iWord))
1655            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1656          else
1657            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1658        }
1659      }
1660    }
1661    way_=1;       // Swap direction
1662  } else {
1663#ifdef FULL_PRINT
1664    printf("Up Fix ");
1665#endif
1666    for (iWord=0;iWord<numberWords;iWord++) {
1667      int i;
1668      for (i=0;i<32;i++) {
1669        unsigned int k = 1<<i;
1670        if ((downMask_[iWord]&k)!=0) {
1671          int iColumn = which[i+32*iWord];
1672#ifdef FULL_PRINT
1673          printf("%d ",i+32*iWord);
1674#endif
1675          // fix weak way
1676          if (clique_->type(i+32*iWord))
1677            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1678          else
1679            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1680        }
1681      }
1682    }
1683    way_=-1;      // Swap direction
1684  }
1685#ifdef FULL_PRINT
1686  printf("\n");
1687#endif
1688  return 0.0;
1689}
1690void
1691CbcLongCliqueBranchingObject::print()
1692{
1693  int iWord;
1694  int numberMembers = clique_->numberMembers();
1695  const int * which = clique_->members();
1696  const int * integerVariables = model_->integerVariable();
1697  int numberWords=(numberMembers+31)>>5;
1698  // *** for way - up means fix all those in down section
1699  if (way_<0) {
1700    printf("Clique - Down Fix ");
1701    for (iWord=0;iWord<numberWords;iWord++) {
1702      int i;
1703      for (i=0;i<32;i++) {
1704        unsigned int k = 1<<i;
1705        if ((upMask_[iWord]&k)!=0) {
1706          int iColumn = which[i+32*iWord];
1707          printf("%d ",integerVariables[iColumn]);
1708        }
1709      }
1710    }
1711  } else {
1712    printf("Clique - Up Fix ");
1713    for (iWord=0;iWord<numberWords;iWord++) {
1714      int i;
1715      for (i=0;i<32;i++) {
1716        unsigned int k = 1<<i;
1717        if ((downMask_[iWord]&k)!=0) {
1718          int iColumn = which[i+32*iWord];
1719          printf("%d ",integerVariables[iColumn]);
1720        }
1721      }
1722    }
1723  }
1724  printf("\n");
1725}
1726// Default Constructor
1727CbcSOSBranchingObject::CbcSOSBranchingObject()
1728  :CbcBranchingObject()
1729{
1730  set_ = NULL;
1731  separator_=0.0;
1732}
1733
1734// Useful constructor
1735CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
1736                                              const CbcSOS * set,
1737                                              int way ,
1738                                              double separator)
1739  :CbcBranchingObject(model,set->id(),way,0.5)
1740{
1741  set_ = set;
1742  separator_ = separator;
1743}
1744
1745// Copy constructor
1746CbcSOSBranchingObject::CbcSOSBranchingObject ( const CbcSOSBranchingObject & rhs) :CbcBranchingObject(rhs)
1747{
1748  set_=rhs.set_;
1749  separator_ = rhs.separator_;
1750}
1751
1752// Assignment operator
1753CbcSOSBranchingObject & 
1754CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject& rhs)
1755{
1756  if (this != &rhs) {
1757    CbcBranchingObject::operator=(rhs);
1758    set_=rhs.set_;
1759    separator_ = rhs.separator_;
1760  }
1761  return *this;
1762}
1763CbcBranchingObject * 
1764CbcSOSBranchingObject::clone() const
1765{ 
1766  return (new CbcSOSBranchingObject(*this));
1767}
1768
1769
1770// Destructor
1771CbcSOSBranchingObject::~CbcSOSBranchingObject ()
1772{
1773}
1774double
1775CbcSOSBranchingObject::branch()
1776{
1777  decrementNumberBranchesLeft();
1778  int numberMembers = set_->numberMembers();
1779  const int * which = set_->members();
1780  const double * weights = set_->weights();
1781  OsiSolverInterface * solver = model_->solver();
1782  //const double * lower = solver->getColLower();
1783  //const double * upper = solver->getColUpper();
1784  // *** for way - up means fix all those in down section
1785  if (way_<0) {
1786    int i;
1787    for ( i=0;i<numberMembers;i++) {
1788      if (weights[i] > separator_)
1789        break;
1790    }
1791    assert (i<numberMembers);
1792    for (;i<numberMembers;i++) 
1793      solver->setColUpper(which[i],0.0);
1794    way_=1;       // Swap direction
1795  } else {
1796    int i;
1797    for ( i=0;i<numberMembers;i++) {
1798      if (weights[i] >= separator_)
1799        break;
1800      else
1801        solver->setColUpper(which[i],0.0);
1802    }
1803    assert (i<numberMembers);
1804    way_=-1;      // Swap direction
1805  }
1806  return 0.0;
1807}
1808// Print what would happen 
1809void
1810CbcSOSBranchingObject::print()
1811{
1812  int numberMembers = set_->numberMembers();
1813  const int * which = set_->members();
1814  const double * weights = set_->weights();
1815  OsiSolverInterface * solver = model_->solver();
1816  //const double * lower = solver->getColLower();
1817  const double * upper = solver->getColUpper();
1818  int first=numberMembers;
1819  int last=-1;
1820  int numberFixed=0;
1821  int numberOther=0;
1822  int i;
1823  for ( i=0;i<numberMembers;i++) {
1824    double bound = upper[which[i]];
1825    if (bound) {
1826      first = CoinMin(first,i);
1827      last = CoinMax(last,i);
1828    }
1829  }
1830  // *** for way - up means fix all those in down section
1831  if (way_<0) {
1832    printf("SOS Down");
1833    for ( i=0;i<numberMembers;i++) {
1834      double bound = upper[which[i]];
1835      if (weights[i] > separator_)
1836        break;
1837      else if (bound)
1838        numberOther++;
1839    }
1840    assert (i<numberMembers);
1841    for (;i<numberMembers;i++) {
1842      double bound = upper[which[i]];
1843      if (bound)
1844        numberFixed++;
1845    }
1846  } else {
1847    printf("SOS Up");
1848    for ( i=0;i<numberMembers;i++) {
1849      double bound = upper[which[i]];
1850      if (weights[i] >= separator_)
1851        break;
1852      else if (bound)
1853        numberFixed++;
1854    }
1855    assert (i<numberMembers);
1856    for (;i<numberMembers;i++) {
1857      double bound = upper[which[i]];
1858      if (bound)
1859        numberOther++;
1860    }
1861  }
1862  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
1863         separator_,which[first],weights[first],which[last],weights[last],numberFixed,numberOther);
1864}
1865 
1866// Default Constructor
1867CbcBranchDefaultDecision::CbcBranchDefaultDecision()
1868  :CbcBranchDecision()
1869{
1870  bestCriterion_ = 0.0;
1871  bestChangeUp_ = 0.0;
1872  bestNumberUp_ = 0;
1873  bestChangeDown_ = 0.0;
1874  bestObject_ = NULL;
1875  bestNumberDown_ = 0;
1876}
1877
1878// Copy constructor
1879CbcBranchDefaultDecision::CbcBranchDefaultDecision (
1880                                    const CbcBranchDefaultDecision & rhs)
1881  :CbcBranchDecision(rhs)
1882{
1883  bestCriterion_ = rhs.bestCriterion_;
1884  bestChangeUp_ = rhs.bestChangeUp_;
1885  bestNumberUp_ = rhs.bestNumberUp_;
1886  bestChangeDown_ = rhs.bestChangeDown_;
1887  bestNumberDown_ = rhs.bestNumberDown_;
1888  bestObject_ = rhs.bestObject_;
1889  model_ = rhs.model_;
1890}
1891
1892CbcBranchDefaultDecision::~CbcBranchDefaultDecision()
1893{
1894}
1895
1896// Clone
1897CbcBranchDecision * 
1898CbcBranchDefaultDecision::clone() const
1899{
1900  return new CbcBranchDefaultDecision(*this);
1901}
1902
1903// Initialize i.e. before start of choosing at a node
1904void 
1905CbcBranchDefaultDecision::initialize(CbcModel * model)
1906{
1907  bestCriterion_ = 0.0;
1908  bestChangeUp_ = 0.0;
1909  bestNumberUp_ = 0;
1910  bestChangeDown_ = 0.0;
1911  bestNumberDown_ = 0;
1912  bestObject_ = NULL;
1913  model_ = model;
1914}
1915
1916
1917/*
1918  Simple default decision algorithm. Compare based on infeasibility (numInfUp,
1919  numInfDn) until a solution is found by search, then switch to change in
1920  objective (changeUp, changeDn). Note that bestSoFar is remembered in
1921  bestObject_, so the parameter bestSoFar is unused.
1922*/
1923
1924int
1925CbcBranchDefaultDecision::betterBranch(CbcBranchingObject * thisOne,
1926                            CbcBranchingObject * bestSoFar,
1927                            double changeUp, int numInfUp,
1928                            double changeDn, int numInfDn)
1929{
1930  bool beforeSolution = cbcModel()->getSolutionCount()==
1931    cbcModel()->getNumberHeuristicSolutions();;
1932  int betterWay=0;
1933  if (beforeSolution) {
1934    if (!bestObject_) {
1935      bestNumberUp_=INT_MAX;
1936      bestNumberDown_=INT_MAX;
1937    }
1938    // before solution - choose smallest number
1939    // could add in depth as well
1940    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1941    if (numInfUp<numInfDn) {
1942      if (numInfUp<bestNumber) {
1943        betterWay = 1;
1944      } else if (numInfUp==bestNumber) {
1945        if (changeUp<bestCriterion_)
1946          betterWay=1;
1947      }
1948    } else if (numInfUp>numInfDn) {
1949      if (numInfDn<bestNumber) {
1950        betterWay = -1;
1951      } else if (numInfDn==bestNumber) {
1952        if (changeDn<bestCriterion_)
1953          betterWay=-1;
1954      }
1955    } else {
1956      // up and down have same number
1957      bool better=false;
1958      if (numInfUp<bestNumber) {
1959        better=true;
1960      } else if (numInfUp==bestNumber) {
1961        if (min(changeUp,changeDn)<bestCriterion_)
1962          better=true;;
1963      }
1964      if (better) {
1965        // see which way
1966        if (changeUp<=changeDn)
1967          betterWay=1;
1968        else
1969          betterWay=-1;
1970      }
1971    }
1972  } else {
1973    if (!bestObject_) {
1974      bestCriterion_=-1.0;
1975    }
1976    // got a solution
1977    if (changeUp<=changeDn) {
1978      if (changeUp>bestCriterion_)
1979        betterWay=1;
1980    } else {
1981      if (changeDn>bestCriterion_)
1982        betterWay=-1;
1983    }
1984  }
1985  if (betterWay) {
1986    bestCriterion_ = CoinMin(changeUp,changeDn);
1987    bestChangeUp_ = changeUp;
1988    bestNumberUp_ = numInfUp;
1989    bestChangeDown_ = changeDn;
1990    bestNumberDown_ = numInfDn;
1991    bestObject_=thisOne;
1992    // See if user is overriding way
1993    if (thisOne->object()&&thisOne->object()->preferredWay())
1994      betterWay = thisOne->object()->preferredWay();
1995  }
1996  return betterWay;
1997}
1998/* Sets or gets best criterion so far */
1999void 
2000CbcBranchDefaultDecision::setBestCriterion(double value)
2001{ 
2002  bestCriterion_ = value;
2003}
2004double 
2005CbcBranchDefaultDecision::getBestCriterion() const
2006{ 
2007  return bestCriterion_;
2008}
2009
2010/* Compare N branching objects. Return index of best
2011   and sets way of branching in chosen object.
2012   
2013   This routine is used only after strong branching.
2014*/
2015
2016int
2017CbcBranchDefaultDecision::bestBranch (CbcBranchingObject ** objects, int numberObjects,
2018                                   int numberUnsatisfied,
2019                                   double * changeUp, int * numberInfeasibilitiesUp,
2020                                   double * changeDown, int * numberInfeasibilitiesDown,
2021                                   double objectiveValue) 
2022{
2023
2024  int bestWay=0;
2025  int whichObject = -1;
2026  if (numberObjects) {
2027    CbcModel * model = cbcModel();
2028    // at continuous
2029    //double continuousObjective = model->getContinuousObjective();
2030    //int continuousInfeasibilities = model->getContinuousInfeasibilities();
2031   
2032    // average cost to get rid of infeasibility
2033    //double averageCostPerInfeasibility =
2034    //(objectiveValue-continuousObjective)/
2035    //(double) (abs(continuousInfeasibilities-numberUnsatisfied)+1);
2036    /* beforeSolution is :
2037       0 - before any solution
2038       n - n heuristic solutions but no branched one
2039       -1 - branched solution found
2040    */
2041    int numberSolutions = model->getSolutionCount();
2042    double cutoff = model->getCutoff();
2043    int method=0;
2044    int i;
2045    if (numberSolutions) {
2046      int numberHeuristic = model->getNumberHeuristicSolutions();
2047      if (numberHeuristic<numberSolutions) {
2048        method = 1;
2049      } else {
2050        method = 2;
2051        // look further
2052        for ( i = 0 ; i < numberObjects ; i++) {
2053          int numberNext = numberInfeasibilitiesUp[i];
2054         
2055          if (numberNext<numberUnsatisfied) {
2056            int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
2057            double perUnsatisfied = changeUp[i]/(double) numberUp;
2058            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2059            if (estimatedObjective<cutoff) 
2060              method=3;
2061          }
2062          numberNext = numberInfeasibilitiesDown[i];
2063          if (numberNext<numberUnsatisfied) {
2064            int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
2065            double perUnsatisfied = changeDown[i]/(double) numberDown;
2066            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2067            if (estimatedObjective<cutoff) 
2068              method=3;
2069          }
2070        }
2071      }
2072      method=2;
2073    } else {
2074      method = 0;
2075    }
2076    // Uncomment next to force method 4
2077    //method=4;
2078    /* Methods :
2079       0 - fewest infeasibilities
2080       1 - largest min change in objective
2081       2 - as 1 but use sum of changes if min close
2082       3 - predicted best solution
2083       4 - take cheapest up branch if infeasibilities same
2084    */
2085    int bestNumber=INT_MAX;
2086    double bestCriterion=-1.0e50;
2087    double alternativeCriterion = -1.0;
2088    double bestEstimate = 1.0e100;
2089    switch (method) {
2090    case 0:
2091      // could add in depth as well
2092      for ( i = 0 ; i < numberObjects ; i++) {
2093        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
2094        if (thisNumber<=bestNumber) {
2095          int betterWay=0;
2096          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
2097            if (numberInfeasibilitiesUp[i]<bestNumber) {
2098              betterWay = 1;
2099            } else {
2100              if (changeUp[i]<bestCriterion)
2101                betterWay=1;
2102            }
2103          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
2104            if (numberInfeasibilitiesDown[i]<bestNumber) {
2105              betterWay = -1;
2106            } else {
2107              if (changeDown[i]<bestCriterion)
2108                betterWay=-1;
2109            }
2110          } else {
2111            // up and down have same number
2112            bool better=false;
2113            if (numberInfeasibilitiesUp[i]<bestNumber) {
2114              better=true;
2115            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2116              if (min(changeUp[i],changeDown[i])<bestCriterion)
2117                better=true;;
2118            }
2119            if (better) {
2120              // see which way
2121              if (changeUp[i]<=changeDown[i])
2122                betterWay=1;
2123              else
2124                betterWay=-1;
2125            }
2126          }
2127          if (betterWay) {
2128            bestCriterion = min(changeUp[i],changeDown[i]);
2129            bestNumber = thisNumber;
2130            whichObject = i;
2131            bestWay = betterWay;
2132          }
2133        }
2134      }
2135      break;
2136    case 1:
2137      for ( i = 0 ; i < numberObjects ; i++) {
2138        int betterWay=0;
2139        if (changeUp[i]<=changeDown[i]) {
2140          if (changeUp[i]>bestCriterion)
2141            betterWay=1;
2142        } else {
2143          if (changeDown[i]>bestCriterion)
2144            betterWay=-1;
2145        }
2146        if (betterWay) {
2147          bestCriterion = min(changeUp[i],changeDown[i]);
2148          whichObject = i;
2149          bestWay = betterWay;
2150        }
2151      }
2152      break;
2153    case 2:
2154      for ( i = 0 ; i < numberObjects ; i++) {
2155        double change = min(changeUp[i],changeDown[i]);
2156        double sum = changeUp[i] + changeDown[i];
2157        bool take=false;
2158        if (change>1.1*bestCriterion) 
2159          take=true;
2160        else if (change>0.9*bestCriterion&&sum+change>bestCriterion+alternativeCriterion) 
2161          take=true;
2162        if (take) {
2163          if (changeUp[i]<=changeDown[i]) {
2164            if (changeUp[i]>bestCriterion)
2165              bestWay=1;
2166          } else {
2167            if (changeDown[i]>bestCriterion)
2168              bestWay=-1;
2169          }
2170          bestCriterion = change;
2171          alternativeCriterion = sum;
2172          whichObject = i;
2173        }
2174      }
2175      break;
2176    case 3:
2177      for ( i = 0 ; i < numberObjects ; i++) {
2178        int numberNext = numberInfeasibilitiesUp[i];
2179       
2180        if (numberNext<numberUnsatisfied) {
2181          int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
2182          double perUnsatisfied = changeUp[i]/(double) numberUp;
2183          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2184          if (estimatedObjective<bestEstimate) {
2185            bestEstimate = estimatedObjective;
2186            bestWay=1;
2187            whichObject=i;
2188          }
2189        }
2190        numberNext = numberInfeasibilitiesDown[i];
2191        if (numberNext<numberUnsatisfied) {
2192          int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
2193          double perUnsatisfied = changeDown[i]/(double) numberDown;
2194          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2195          if (estimatedObjective<bestEstimate) {
2196            bestEstimate = estimatedObjective;
2197            bestWay=-1;
2198            whichObject=i;
2199          }
2200        }
2201      }
2202      break;
2203    case 4:
2204      // if number infeas same then cheapest up
2205      // first get best number or when going down
2206      // now choose smallest change up amongst equal number infeas
2207      for ( i = 0 ; i < numberObjects ; i++) {
2208        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
2209        if (thisNumber<=bestNumber) {
2210          int betterWay=0;
2211          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
2212            if (numberInfeasibilitiesUp[i]<bestNumber) {
2213              betterWay = 1;
2214            } else {
2215              if (changeUp[i]<bestCriterion)
2216                betterWay=1;
2217            }
2218          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
2219            if (numberInfeasibilitiesDown[i]<bestNumber) {
2220              betterWay = -1;
2221            } else {
2222              if (changeDown[i]<bestCriterion)
2223                betterWay=-1;
2224            }
2225          } else {
2226            // up and down have same number
2227            bool better=false;
2228            if (numberInfeasibilitiesUp[i]<bestNumber) {
2229              better=true;
2230            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2231              if (min(changeUp[i],changeDown[i])<bestCriterion)
2232                better=true;;
2233            }
2234            if (better) {
2235              // see which way
2236              if (changeUp[i]<=changeDown[i])
2237                betterWay=1;
2238              else
2239                betterWay=-1;
2240            }
2241          }
2242          if (betterWay) {
2243            bestCriterion = min(changeUp[i],changeDown[i]);
2244            bestNumber = thisNumber;
2245            whichObject = i;
2246            bestWay = betterWay;
2247          }
2248        }
2249      }
2250      bestCriterion=1.0e50;
2251      for ( i = 0 ; i < numberObjects ; i++) {
2252        int thisNumber = numberInfeasibilitiesUp[i];
2253        if (thisNumber==bestNumber&&changeUp) {
2254          if (changeUp[i]<bestCriterion) {
2255            bestCriterion = changeUp[i];
2256            whichObject = i;
2257            bestWay = 1;
2258          }
2259        }
2260      }
2261      break;
2262    }
2263    // set way in best
2264    if (whichObject>=0) {
2265      CbcBranchingObject * bestObject = objects[whichObject];
2266      if (bestObject->object()&&bestObject->object()->preferredWay()) 
2267        bestWay = bestObject->object()->preferredWay();
2268      bestObject->way(bestWay);
2269    } else {
2270      printf("debug\n");
2271    }
2272  }
2273  return whichObject;
2274}
2275
2276// Default Constructor
2277CbcFollowOn::CbcFollowOn ()
2278  : CbcObject(),
2279    rhs_(NULL)
2280{
2281}
2282
2283// Useful constructor
2284CbcFollowOn::CbcFollowOn (CbcModel * model)
2285  : CbcObject(model)
2286{
2287  assert (model);
2288  OsiSolverInterface * solver = model_->solver();
2289  matrix_ = *solver->getMatrixByCol();
2290  matrix_.removeGaps();
2291  matrix_.setExtraGap(0.0);
2292  matrixByRow_ = *solver->getMatrixByRow();
2293  int numberRows = matrix_.getNumRows();
2294 
2295  rhs_ = new int[numberRows];
2296  int i;
2297  const double * rowLower = solver->getRowLower();
2298  const double * rowUpper = solver->getRowUpper();
2299  // Row copy
2300  const double * elementByRow = matrixByRow_.getElements();
2301  const int * column = matrixByRow_.getIndices();
2302  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2303  const int * rowLength = matrixByRow_.getVectorLengths();
2304  for (i=0;i<numberRows;i++) {
2305    rhs_[i]=0;
2306    double value = rowLower[i];
2307    if (value==rowUpper[i]) {
2308      if (floor(value)==value&&value>=1.0&&value<10.0) {
2309        // check elements
2310        bool good=true;
2311        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2312          int iColumn = column[j];
2313          if (!solver->isBinary(iColumn))
2314            good=false;
2315          double elValue = elementByRow[j];
2316          if (floor(elValue)!=elValue||value<1.0)
2317            good=false;
2318        }
2319        if (good)
2320          rhs_[i]=(int) value;
2321      }
2322    }
2323  }
2324}
2325
2326// Copy constructor
2327CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
2328  :CbcObject(rhs),
2329   matrix_(rhs.matrix_),
2330   matrixByRow_(rhs.matrixByRow_)
2331{
2332  int numberRows = matrix_.getNumRows();
2333  rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2334}
2335
2336// Clone
2337CbcObject *
2338CbcFollowOn::clone() const
2339{
2340  return new CbcFollowOn(*this);
2341}
2342
2343// Assignment operator
2344CbcFollowOn & 
2345CbcFollowOn::operator=( const CbcFollowOn& rhs)
2346{
2347  if (this!=&rhs) {
2348    CbcObject::operator=(rhs);
2349    delete [] rhs_;
2350    matrix_ = rhs.matrix_;
2351    matrixByRow_ = rhs.matrixByRow_;
2352    int numberRows = matrix_.getNumRows();
2353    rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2354  }
2355  return *this;
2356}
2357
2358// Destructor
2359CbcFollowOn::~CbcFollowOn ()
2360{
2361  delete [] rhs_;
2362}
2363// As some computation is needed in more than one place - returns row
2364int 
2365CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
2366{
2367  int whichRow=-1;
2368  otherRow=-1;
2369  int numberRows = matrix_.getNumRows();
2370 
2371  int i;
2372  // For sorting
2373  int * sort = new int [numberRows];
2374  int * isort = new int [numberRows];
2375  // Column copy
2376  //const double * element = matrix_.getElements();
2377  const int * row = matrix_.getIndices();
2378  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2379  const int * columnLength = matrix_.getVectorLengths();
2380  // Row copy
2381  const double * elementByRow = matrixByRow_.getElements();
2382  const int * column = matrixByRow_.getIndices();
2383  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2384  const int * rowLength = matrixByRow_.getVectorLengths();
2385  OsiSolverInterface * solver = model_->solver();
2386  const double * columnLower = solver->getColLower();
2387  const double * columnUpper = solver->getColUpper();
2388  const double * solution = solver->getColSolution();
2389  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2390  int nSort=0;
2391  for (i=0;i<numberRows;i++) {
2392    if (rhs_[i]) {
2393      // check elements
2394      double smallest=1.0e10;
2395      double largest=0.0;
2396      int rhsValue=rhs_[i];
2397      int number1=0;
2398      int numberUnsatisfied=0;
2399      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2400        int iColumn = column[j];
2401        double value = elementByRow[j];
2402        double solValue = solution[iColumn];
2403        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2404          smallest = CoinMin(smallest,value);
2405          largest = CoinMax(largest,value);
2406          if (value==1.0)
2407            number1++;
2408          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
2409            numberUnsatisfied++;
2410        } else {
2411          rhsValue -= (int)(value*floor(solValue+0.5));
2412        }
2413      }
2414      if (numberUnsatisfied>1) {
2415        if (smallest<largest) {
2416          // probably no good but check a few things
2417          assert (largest<=rhsValue);
2418          if (number1==1&&largest==rhsValue)
2419            printf("could fix\n");
2420        } else if (largest==rhsValue) {
2421          sort[nSort]=i;
2422          isort[nSort++]=-numberUnsatisfied;
2423        }
2424      }
2425    }
2426  }
2427  if (nSort>1) {
2428    CoinSort_2(isort,isort+nSort,sort);
2429    CoinZeroN(isort,numberRows);
2430    double * other = new double[numberRows];
2431    CoinZeroN(other,numberRows);
2432    int * which = new int[numberRows];
2433    //#define COUNT
2434#ifndef COUNT
2435    bool beforeSolution = model_->getSolutionCount()==0;
2436#endif
2437    for (int k=0;k<nSort-1;k++) {
2438      i=sort[k];
2439      int numberUnsatisfied = 0;
2440      int n=0;
2441      int j;
2442      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2443        int iColumn = column[j];
2444        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2445          double solValue = solution[iColumn]-columnLower[iColumn];
2446          if (solValue<1.0-integerTolerance&&solValue>integerTolerance) {
2447            numberUnsatisfied++;
2448            for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2449              int iRow = row[jj];
2450              if (rhs_[iRow]) {
2451                other[iRow]+=solValue;
2452                if (isort[iRow]) {
2453                  isort[iRow]++;
2454                } else {
2455                  isort[iRow]=1;
2456                  which[n++]=iRow;
2457                }
2458              }
2459            }
2460          }
2461        }
2462      }
2463      double total=0.0;
2464      // Take out row
2465      double sumThis=other[i];
2466      other[i]=0.0;
2467      assert (numberUnsatisfied==isort[i]);
2468      // find one nearest half if solution, one if before solution
2469      int iBest=-1;
2470      double dtarget=0.5*total;
2471#ifdef COUNT
2472      int target = (numberUnsatisfied+1)>>1;
2473      int best=numberUnsatisfied;
2474#else
2475      double best;
2476      if (beforeSolution)
2477        best=dtarget;
2478      else
2479        best=1.0e30;
2480#endif
2481      for (j=0;j<n;j++) {
2482        int iRow = which[j];
2483        double dvalue=other[iRow];
2484        other[iRow]=0.0;
2485#ifdef COUNT
2486        int value = isort[iRow];
2487#endif
2488        isort[iRow]=0;
2489        if (fabs(dvalue)<1.0e-8||fabs(sumThis-dvalue)<1.0e-8)
2490          continue;
2491        if (dvalue<integerTolerance||dvalue>1.0-integerTolerance)
2492          continue;
2493#ifdef COUNT
2494        if (abs(value-target)<best&&value!=numberUnsatisfied) {
2495          best=abs(value-target);
2496          iBest=iRow;
2497          if (dvalue<dtarget)
2498            preferredWay=1;
2499          else
2500            preferredWay=-1;
2501        }
2502#else
2503        if (beforeSolution) {
2504          if (fabs(dvalue-dtarget)>best) {
2505            best = fabs(dvalue-dtarget);
2506            iBest=iRow;
2507            if (dvalue<dtarget)
2508              preferredWay=1;
2509            else
2510              preferredWay=-1;
2511          }
2512        } else {
2513          if (fabs(dvalue-dtarget)<best) {
2514            best = fabs(dvalue-dtarget);
2515            iBest=iRow;
2516            if (dvalue<dtarget)
2517              preferredWay=1;
2518            else
2519              preferredWay=-1;
2520          }
2521        }
2522#endif
2523      }
2524      if (iBest>=0) {
2525        whichRow=i;
2526        otherRow=iBest;
2527        break;
2528      }
2529    }
2530    delete [] which;
2531    delete [] other;
2532  }
2533  delete [] sort;
2534  delete [] isort;
2535  return whichRow;
2536}
2537
2538// Infeasibility - large is 0.5
2539double 
2540CbcFollowOn::infeasibility(int & preferredWay) const
2541{
2542  int otherRow=0;
2543  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2544  if (whichRow<0)
2545    return 0.0;
2546  else
2547  return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
2548}
2549
2550// This looks at solution and sets bounds to contain solution
2551void 
2552CbcFollowOn::feasibleRegion()
2553{
2554}
2555
2556
2557// Creates a branching object
2558CbcBranchingObject * 
2559CbcFollowOn::createBranch(int way) 
2560{
2561  int otherRow=0;
2562  int preferredWay;
2563  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2564  assert(way==preferredWay);
2565  assert (whichRow>=0);
2566  int numberColumns = matrix_.getNumCols();
2567 
2568  // Column copy
2569  //const double * element = matrix_.getElements();
2570  const int * row = matrix_.getIndices();
2571  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2572  const int * columnLength = matrix_.getVectorLengths();
2573  // Row copy
2574  //const double * elementByRow = matrixByRow_.getElements();
2575  const int * column = matrixByRow_.getIndices();
2576  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2577  const int * rowLength = matrixByRow_.getVectorLengths();
2578  OsiSolverInterface * solver = model_->solver();
2579  const double * columnLower = solver->getColLower();
2580  const double * columnUpper = solver->getColUpper();
2581  //const double * solution = solver->getColSolution();
2582  int nUp=0;
2583  int nDown=0;
2584  int * upList = new int[numberColumns];
2585  int * downList = new int[numberColumns];
2586  int j;
2587  for (j=rowStart[whichRow];j<rowStart[whichRow]+rowLength[whichRow];j++) {
2588    int iColumn = column[j];
2589    if (columnLower[iColumn]!=columnUpper[iColumn]) {
2590      bool up=true;
2591      for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2592        int iRow = row[jj];
2593        if (iRow==otherRow) {
2594          up=false;
2595          break;
2596        }
2597      }
2598      if (up)
2599        upList[nUp++]=iColumn;
2600      else
2601        downList[nDown++]=iColumn;
2602    }
2603  }
2604  //printf("way %d\n",way);
2605  // create object
2606  //printf("would fix %d down and %d up\n",nDown,nUp);
2607  CbcBranchingObject * branch
2608     = new CbcFixingBranchingObject(model_,way,
2609                                         nDown,downList,nUp,upList);
2610  delete [] upList;
2611  delete [] downList;
2612  return branch;
2613}
2614// Default Constructor
2615CbcFixingBranchingObject::CbcFixingBranchingObject()
2616  :CbcBranchingObject()
2617{
2618  numberDown_=0;
2619  numberUp_=0;
2620  downList_=NULL;
2621  upList_=NULL;
2622}
2623
2624// Useful constructor
2625CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
2626                                                    int way ,
2627                                                    int numberOnDownSide, const int * down,
2628                                                    int numberOnUpSide, const int * up)
2629  :CbcBranchingObject(model,0,way,0.5)
2630{
2631  numberDown_=numberOnDownSide;
2632  numberUp_=numberOnUpSide;
2633  downList_ = CoinCopyOfArray(down,numberDown_);
2634  upList_ = CoinCopyOfArray(up,numberUp_);
2635}
2636
2637// Copy constructor
2638CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) :CbcBranchingObject(rhs)
2639{
2640  numberDown_=rhs.numberDown_;
2641  numberUp_=rhs.numberUp_;
2642  downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2643  upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2644}
2645
2646// Assignment operator
2647CbcFixingBranchingObject & 
2648CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject& rhs)
2649{
2650  if (this != &rhs) {
2651    CbcBranchingObject::operator=(rhs);
2652    delete [] downList_;
2653    delete [] upList_;
2654    numberDown_=rhs.numberDown_;
2655    numberUp_=rhs.numberUp_;
2656    downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2657    upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2658  }
2659  return *this;
2660}
2661CbcBranchingObject * 
2662CbcFixingBranchingObject::clone() const
2663{ 
2664  return (new CbcFixingBranchingObject(*this));
2665}
2666
2667
2668// Destructor
2669CbcFixingBranchingObject::~CbcFixingBranchingObject ()
2670{
2671  delete [] downList_;
2672  delete [] upList_;
2673}
2674double
2675CbcFixingBranchingObject::branch()
2676{
2677  decrementNumberBranchesLeft();
2678  OsiSolverInterface * solver = model_->solver();
2679  const double * columnLower = solver->getColLower();
2680  int i;
2681  // *** for way - up means fix all those in up section
2682  if (way_<0) {
2683#ifdef FULL_PRINT
2684    printf("Down Fix ");
2685#endif
2686    //printf("Down Fix %d\n",numberDown_);
2687    for (i=0;i<numberDown_;i++) {
2688      int iColumn = downList_[i];
2689      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2690#ifdef FULL_PRINT
2691      printf("Setting bound on %d to lower bound\n",iColumn);
2692#endif
2693    }
2694    way_=1;       // Swap direction
2695  } else {
2696#ifdef FULL_PRINT
2697    printf("Up Fix ");
2698#endif
2699    //printf("Up Fix %d\n",numberUp_);
2700    for (i=0;i<numberUp_;i++) {
2701      int iColumn = upList_[i];
2702      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2703#ifdef FULL_PRINT
2704      printf("Setting bound on %d to lower bound\n",iColumn);
2705#endif
2706    }
2707    way_=-1;      // Swap direction
2708  }
2709#ifdef FULL_PRINT
2710  printf("\n");
2711#endif
2712  return 0.0;
2713}
2714void
2715CbcFixingBranchingObject::print()
2716{
2717  int i;
2718  // *** for way - up means fix all those in up section
2719  if (way_<0) {
2720    printf("Down Fix ");
2721    for (i=0;i<numberDown_;i++) {
2722      int iColumn = downList_[i];
2723      printf("%d ",iColumn);
2724    }
2725  } else {
2726    printf("Up Fix ");
2727    for (i=0;i<numberUp_;i++) {
2728      int iColumn = upList_[i];
2729      printf("%d ",iColumn);
2730    }
2731  }
2732  printf("\n");
2733}
2734// Default Constructor
2735CbcNWay::CbcNWay ()
2736  : CbcObject(),
2737    numberMembers_(0),
2738    members_(NULL),
2739    consequence_(NULL)
2740{
2741}
2742
2743// Useful constructor (which are integer indices)
2744CbcNWay::CbcNWay (CbcModel * model, int numberMembers,
2745                  const int * which, int identifier)
2746  : CbcObject(model)
2747{
2748  id_=identifier;
2749  numberMembers_=numberMembers;
2750  if (numberMembers_) {
2751    members_ = new int[numberMembers_];
2752    memcpy(members_,which,numberMembers_*sizeof(int));
2753  } else {
2754    members_ = NULL;
2755  }
2756  consequence_ = NULL;
2757}
2758
2759// Copy constructor
2760CbcNWay::CbcNWay ( const CbcNWay & rhs)
2761  :CbcObject(rhs)
2762{
2763  numberMembers_ = rhs.numberMembers_;
2764  consequence_ = NULL;
2765  if (numberMembers_) {
2766    members_ = new int[numberMembers_];
2767    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
2768    if (rhs.consequence_) {
2769      consequence_ = new CbcConsequence * [numberMembers_];
2770      for (int i=0;i<numberMembers_;i++) {
2771        if (rhs.consequence_[i])
2772          consequence_[i]= rhs.consequence_[i]->clone();
2773        else
2774          consequence_[i]=NULL;
2775      }
2776    }
2777  } else {
2778    members_ = NULL;
2779  }
2780}
2781
2782// Clone
2783CbcObject *
2784CbcNWay::clone() const
2785{
2786  return new CbcNWay(*this);
2787}
2788
2789// Assignment operator
2790CbcNWay & 
2791CbcNWay::operator=( const CbcNWay& rhs)
2792{
2793  if (this!=&rhs) {
2794    CbcObject::operator=(rhs);
2795    delete [] members_;
2796    numberMembers_ = rhs.numberMembers_;
2797    if (consequence_) {
2798      for (int i=0;i<numberMembers_;i++) 
2799        delete consequence_[i];
2800      delete [] consequence_;
2801      consequence_=NULL;
2802    }
2803    if (numberMembers_) {
2804      members_ = new int[numberMembers_];
2805      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
2806    } else {
2807      members_ = NULL;
2808    }
2809    if (rhs.consequence_) {
2810      consequence_ = new CbcConsequence * [numberMembers_];
2811      for (int i=0;i<numberMembers_;i++) {
2812        if (rhs.consequence_[i])
2813          consequence_[i]= rhs.consequence_[i]->clone();
2814        else
2815          consequence_[i]=NULL;
2816      }
2817    }
2818  }
2819  return *this;
2820}
2821
2822// Destructor
2823CbcNWay::~CbcNWay ()
2824{
2825  delete [] members_;
2826  if (consequence_) {
2827    for (int i=0;i<numberMembers_;i++) 
2828      delete consequence_[i];
2829    delete [] consequence_;
2830  }
2831}
2832// Set up a consequence for a single member
2833void 
2834CbcNWay::setConsequence(int iColumn, const CbcConsequence & consequence)
2835{
2836  if (!consequence_) {
2837    consequence_ = new CbcConsequence * [numberMembers_];
2838    for (int i=0;i<numberMembers_;i++) 
2839      consequence_[i]=NULL;
2840  }
2841  for (int i=0;i<numberMembers_;i++) {
2842    if (members_[i]==iColumn) {
2843      consequence_[i]=consequence.clone();
2844      break;
2845    }
2846  }
2847}
2848
2849// Applies a consequence for a single member
2850void 
2851CbcNWay::applyConsequence(int iSequence, int state) const
2852{
2853  assert (state==-9999||state==9999);
2854  if (consequence_) {
2855    CbcConsequence * consequence = consequence_[iSequence];
2856    if (consequence) 
2857      consequence->applyToSolver(model_->solver(),state);
2858  }
2859}
2860 
2861// Infeasibility - large is 0.5
2862double 
2863CbcNWay::infeasibility(int & preferredWay) const
2864{
2865  int numberUnsatis=0;
2866  int j;
2867  OsiSolverInterface * solver = model_->solver();
2868  const double * solution = model_->testSolution();
2869  const double * lower = solver->getColLower();
2870  const double * upper = solver->getColUpper();
2871  double largestValue=0.0;
2872 
2873  double integerTolerance = 
2874    model_->getDblParam(CbcModel::CbcIntegerTolerance);
2875
2876  for (j=0;j<numberMembers_;j++) {
2877    int iColumn = members_[j];
2878    double value = solution[iColumn];
2879    value = CoinMax(value, lower[iColumn]);
2880    value = CoinMin(value, upper[iColumn]);
2881    double distance = CoinMin(value-lower[iColumn],upper[iColumn]-value);
2882    if (distance>integerTolerance) {
2883      numberUnsatis++;
2884      largestValue = CoinMax(distance,largestValue);
2885    }
2886  }
2887  preferredWay=1;
2888  if (numberUnsatis) {
2889    return largestValue;
2890  } else {
2891    return 0.0; // satisfied
2892  }
2893}
2894
2895// This looks at solution and sets bounds to contain solution
2896void 
2897CbcNWay::feasibleRegion()
2898{
2899  int j;
2900  OsiSolverInterface * solver = model_->solver();
2901  const double * solution = model_->testSolution();
2902  const double * lower = solver->getColLower();
2903  const double * upper = solver->getColUpper();
2904  double integerTolerance = 
2905    model_->getDblParam(CbcModel::CbcIntegerTolerance);
2906  for (j=0;j<numberMembers_;j++) {
2907    int iColumn = members_[j];
2908    double value = solution[iColumn];
2909    value = CoinMax(value, lower[iColumn]);
2910    value = CoinMin(value, upper[iColumn]);
2911    if (value>=upper[iColumn]-integerTolerance) {
2912      solver->setColLower(iColumn,upper[iColumn]);
2913    } else {
2914      assert (value<=lower[iColumn]+integerTolerance);
2915      solver->setColUpper(iColumn,lower[iColumn]);
2916    }
2917  }
2918}
2919// Redoes data when sequence numbers change
2920void 
2921CbcNWay::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
2922{
2923  model_=model;
2924  int n2=0;
2925  for (int j=0;j<numberMembers_;j++) {
2926    int iColumn = members_[j];
2927    int i;
2928    for (i=0;i<numberColumns;i++) {
2929      if (originalColumns[i]==iColumn)
2930        break;
2931    }
2932    if (i<numberColumns) {
2933      members_[n2]=i;
2934      consequence_[n2++]=consequence_[j];
2935    } else {
2936      delete consequence_[j];
2937    }
2938  }
2939  if (n2<numberMembers_) {
2940    printf("** NWay number of members reduced from %d to %d!\n",numberMembers_,n2);
2941    numberMembers_=n2;
2942  }
2943}
2944
2945
2946// Creates a branching object
2947CbcBranchingObject * 
2948CbcNWay::createBranch(int way) 
2949{
2950  int numberFree=0;
2951  int j;
2952
2953  OsiSolverInterface * solver = model_->solver();
2954  const double * solution = model_->testSolution();
2955  const double * lower = solver->getColLower();
2956  const double * upper = solver->getColUpper();
2957  int * list = new int[numberMembers_];
2958  double * sort = new double[numberMembers_];
2959
2960  for (j=0;j<numberMembers_;j++) {
2961    int iColumn = members_[j];
2962    double value = solution[iColumn];
2963    value = CoinMax(value, lower[iColumn]);
2964    value = CoinMin(value, upper[iColumn]);
2965    if (upper[iColumn]>lower[iColumn]) {
2966      double distance = upper[iColumn]-value;
2967      list[numberFree]=j;
2968      sort[numberFree++]=distance;
2969    }
2970  }
2971  assert (numberFree);
2972  // sort
2973  CoinSort_2(sort,sort+numberFree,list);
2974  // create object
2975  CbcBranchingObject * branch;
2976  branch = new CbcNWayBranchingObject(model_,this,numberFree,list);
2977  branch->setOriginalObject(this);
2978  delete [] list;
2979  delete [] sort;
2980  return branch;
2981}
2982 
2983// Default Constructor
2984CbcNWayBranchingObject::CbcNWayBranchingObject()
2985  :CbcBranchingObject()
2986{
2987  order_=NULL;
2988  object_=NULL;
2989  numberInSet_=0;
2990  way_=0;
2991}
2992
2993// Useful constructor
2994CbcNWayBranchingObject::CbcNWayBranchingObject (CbcModel * model,
2995                                                const CbcNWay * nway, 
2996                                                int number, const int * order)
2997  :CbcBranchingObject(model,nway->id(),-1,0.5)
2998{
2999  numberBranches_ = number;
3000  order_ = new int [number];
3001  object_=nway;
3002  numberInSet_=number;
3003  memcpy(order_,order,number*sizeof(int));
3004}
3005
3006// Copy constructor
3007CbcNWayBranchingObject::CbcNWayBranchingObject ( const CbcNWayBranchingObject & rhs) :CbcBranchingObject(rhs)
3008{
3009  numberInSet_=rhs.numberInSet_;
3010  object_=rhs.object_;
3011  if (numberInSet_) {
3012    order_ = new int [numberInSet_];
3013    memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
3014  } else {
3015    order_=NULL;
3016  }   
3017}
3018
3019// Assignment operator
3020CbcNWayBranchingObject & 
3021CbcNWayBranchingObject::operator=( const CbcNWayBranchingObject& rhs)
3022{
3023  if (this != &rhs) {
3024    CbcBranchingObject::operator=(rhs);
3025    object_=rhs.object_;
3026    delete [] order_;
3027    numberInSet_=rhs.numberInSet_;
3028    if (numberInSet_) {
3029      order_ = new int [numberInSet_];
3030      memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
3031    } else {
3032      order_=NULL;
3033    }   
3034  }
3035  return *this;
3036}
3037CbcBranchingObject * 
3038CbcNWayBranchingObject::clone() const
3039{ 
3040  return (new CbcNWayBranchingObject(*this));
3041}
3042
3043
3044// Destructor
3045CbcNWayBranchingObject::~CbcNWayBranchingObject ()
3046{
3047  delete [] order_;
3048}
3049double
3050CbcNWayBranchingObject::branch()
3051{
3052  int which = branchIndex_;
3053  branchIndex_++;
3054  assert (numberBranchesLeft()>=0);
3055  if (which==0) {
3056    // first branch so way_ may mean something
3057    assert (way_==-1||way_==1);
3058    if (way_==-1)
3059      which++;
3060  } else if (which==1) {
3061    // second branch so way_ may mean something
3062    assert (way_==-1||way_==1);
3063    if (way_==-1)
3064      which--;
3065    // switch way off
3066    way_=0;
3067  }
3068  const double * lower = model_->solver()->getColLower();
3069  const double * upper = model_->solver()->getColUpper();
3070  const int * members = object_->members();
3071  for (int j=0;j<numberInSet_;j++) {
3072    int iSequence = order_[j];
3073    int iColumn = members[iSequence];
3074    if (j!=which) {
3075      model_->solver()->setColUpper(iColumn,lower[iColumn]);
3076      //model_->solver()->setColLower(iColumn,lower[iColumn]);
3077      assert (lower[iColumn]>-1.0e20);
3078      // apply any consequences
3079      object_->applyConsequence(iSequence,-9999);
3080    } else {
3081      model_->solver()->setColLower(iColumn,upper[iColumn]);
3082      //model_->solver()->setColUpper(iColumn,upper[iColumn]);
3083#ifdef FULL_PRINT
3084      printf("Up Fix %d to %g\n",iColumn,upper[iColumn]);
3085#endif
3086      assert (upper[iColumn]<1.0e20);
3087      // apply any consequences
3088      object_->applyConsequence(iSequence,9999);
3089    }
3090  }
3091  return 0.0;
3092}
3093void
3094CbcNWayBranchingObject::print()
3095{
3096  printf("NWay - Up Fix ");
3097  const int * members = object_->members();
3098  for (int j=0;j<way_;j++) {
3099    int iColumn = members[order_[j]];
3100    printf("%d ",iColumn);
3101  }
3102  printf("\n");
3103}
3104
3105// Default Constructor
3106CbcFixVariable::CbcFixVariable ()
3107  : CbcConsequence(),
3108    numberStates_(0),
3109    states_(NULL),
3110    startLower_(NULL),
3111    startUpper_(NULL),
3112    newBound_(NULL),
3113    variable_(NULL)
3114{
3115}
3116
3117// One useful Constructor
3118CbcFixVariable::CbcFixVariable (int numberStates,const int * states, const int * numberNewLower, 
3119                                const int ** newLowerValue,
3120                                const int ** lowerColumn,
3121                                const int * numberNewUpper, const int ** newUpperValue,
3122                                const int ** upperColumn)
3123  : CbcConsequence(),
3124    states_(NULL),
3125    startLower_(NULL),
3126    startUpper_(NULL),
3127    newBound_(NULL),
3128    variable_(NULL)
3129{
3130  // How much space
3131  numberStates_ = numberStates;
3132  if (numberStates_) {
3133    states_ = new int[numberStates_];
3134    memcpy(states_,states,numberStates_*sizeof(int));
3135    int i;
3136    int n=0;
3137    startLower_ = new int[numberStates_+1];
3138    startUpper_ = new int[numberStates_+1];
3139    startLower_[0]=0;
3140    //count
3141    for (i=0;i<numberStates_;i++) {
3142      n += numberNewLower[i];
3143      startUpper_[i]=n;
3144      n += numberNewUpper[i];
3145      startLower_[i+1]=n;
3146    }
3147    newBound_ = new double [n];
3148    variable_ = new int [n];
3149    n=0;
3150    for (i=0;i<numberStates_;i++) {
3151      int j;
3152      int k;
3153      const int * bound;
3154      const int * variable;
3155      k=numberNewLower[i];
3156      bound = newLowerValue[i];
3157      variable = lowerColumn[i];
3158      for (j=0;j<k;j++) {
3159        newBound_[n]=bound[j];
3160        variable_[n++]=variable[j];
3161      }
3162      k=numberNewUpper[i];
3163      bound = newUpperValue[i];
3164      variable = upperColumn[i];
3165      for (j=0;j<k;j++) {
3166        newBound_[n]=bound[j];
3167        variable_[n++]=variable[j];
3168      }
3169    }
3170  }
3171}
3172
3173// Copy constructor
3174CbcFixVariable::CbcFixVariable ( const CbcFixVariable & rhs)
3175  :CbcConsequence(rhs)
3176{
3177  numberStates_ = rhs.numberStates_;
3178  states_ = NULL;
3179  startLower_ = NULL;
3180  startUpper_ = NULL;
3181  newBound_ = NULL;
3182  variable_ = NULL;
3183  if (numberStates_) {
3184    states_ = CoinCopyOfArray(rhs.states_,numberStates_);
3185    startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
3186    startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
3187    int n=startLower_[numberStates_];
3188    newBound_ = CoinCopyOfArray(rhs.newBound_,n);
3189    variable_ = CoinCopyOfArray(rhs.variable_,n);
3190  }
3191}
3192
3193// Clone
3194CbcConsequence *
3195CbcFixVariable::clone() const
3196{
3197  return new CbcFixVariable(*this);
3198}
3199
3200// Assignment operator
3201CbcFixVariable & 
3202CbcFixVariable::operator=( const CbcFixVariable& rhs)
3203{
3204  if (this!=&rhs) {
3205    CbcConsequence::operator=(rhs);
3206    delete [] states_;
3207    delete [] startLower_;
3208    delete [] startUpper_;
3209    delete [] newBound_;
3210    delete [] variable_;
3211    states_ = NULL;
3212    startLower_ = NULL;
3213    startUpper_ = NULL;
3214    newBound_ = NULL;
3215    variable_ = NULL;
3216    numberStates_ = rhs.numberStates_;
3217    if (numberStates_) {
3218      states_ = CoinCopyOfArray(rhs.states_,numberStates_);
3219      startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
3220      startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
3221      int n=startLower_[numberStates_];
3222      newBound_ = CoinCopyOfArray(rhs.newBound_,n);
3223      variable_ = CoinCopyOfArray(rhs.variable_,n);
3224    }
3225  }
3226  return *this;
3227}
3228
3229// Destructor
3230CbcFixVariable::~CbcFixVariable ()
3231{
3232  delete [] states_;
3233  delete [] startLower_;
3234  delete [] startUpper_;
3235  delete [] newBound_;
3236  delete [] variable_;
3237}
3238// Set up a startLower for a single member
3239void 
3240CbcFixVariable::applyToSolver(OsiSolverInterface * solver, int state) const
3241{
3242  assert (state==-9999||state==9999);
3243  // Find state
3244  int find;
3245  for (find=0;find<numberStates_;find++) 
3246    if (states_[find]==state)
3247      break;
3248  if (find==numberStates_)
3249    return;
3250  int i;
3251  // Set new lower bounds
3252  for (i=startLower_[find];i<startUpper_[find];i++) {
3253    int iColumn = variable_[i];
3254    double value = newBound_[i];
3255    double oldValue = solver->getColLower()[iColumn];
3256    //printf("for %d old lower bound %g, new %g",iColumn,oldValue,value);
3257    solver->setColLower(iColumn,CoinMax(value,oldValue));
3258    //printf(" => %g\n",solver->getColLower()[iColumn]);
3259  }
3260  // Set new upper bounds
3261  for (i=startUpper_[find];i<startLower_[find+1];i++) {
3262    int iColumn = variable_[i];
3263    double value = newBound_[i];
3264    double oldValue = solver->getColUpper()[iColumn];
3265    //printf("for %d old upper bound %g, new %g",iColumn,oldValue,value);
3266    solver->setColUpper(iColumn,CoinMin(value,oldValue));
3267    //printf(" => %g\n",solver->getColUpper()[iColumn]);
3268  }
3269}
3270
3271// Default Constructor
3272CbcDummyBranchingObject::CbcDummyBranchingObject(CbcModel * model)
3273  :CbcBranchingObject(model,0,0,0.5)
3274{
3275  setNumberBranchesLeft(1);
3276}
3277
3278
3279// Copy constructor
3280CbcDummyBranchingObject::CbcDummyBranchingObject ( const CbcDummyBranchingObject & rhs) :CbcBranchingObject(rhs)
3281{
3282}
3283
3284// Assignment operator
3285CbcDummyBranchingObject & 
3286CbcDummyBranchingObject::operator=( const CbcDummyBranchingObject& rhs)
3287{
3288  if (this != &rhs) {
3289    CbcBranchingObject::operator=(rhs);
3290  }
3291  return *this;
3292}
3293CbcBranchingObject * 
3294CbcDummyBranchingObject::clone() const
3295{ 
3296  return (new CbcDummyBranchingObject(*this));
3297}
3298
3299
3300// Destructor
3301CbcDummyBranchingObject::~CbcDummyBranchingObject ()
3302{
3303}
3304
3305/*
3306  Perform a dummy branch
3307*/
3308double
3309CbcDummyBranchingObject::branch()
3310{
3311  decrementNumberBranchesLeft();
3312  return 0.0;
3313}
3314// Print what would happen 
3315void
3316CbcDummyBranchingObject::print()
3317{
3318  printf("Dummy branch\n");
3319}
Note: See TracBrowser for help on using the repository browser.