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

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

for osi methods

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 86.0 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#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());
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  if (way_<0) {
988#ifdef CBC_DEBUG
989  { double olb,oub ;
990    olb = model_->solver()->getColLower()[iColumn] ;
991    oub = model_->solver()->getColUpper()[iColumn] ;
992    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
993           iColumn,olb,oub,down_[0],down_[1]) ; }
994#endif
995    model_->solver()->setColLower(iColumn,down_[0]);
996    model_->solver()->setColUpper(iColumn,down_[1]);
997    way_=1;
998  } else {
999#ifdef CBC_DEBUG
1000  { double olb,oub ;
1001    olb = model_->solver()->getColLower()[iColumn] ;
1002    oub = model_->solver()->getColUpper()[iColumn] ;
1003    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
1004           iColumn,olb,oub,up_[0],up_[1]) ; }
1005#endif
1006    model_->solver()->setColLower(iColumn,up_[0]);
1007    model_->solver()->setColUpper(iColumn,up_[1]);
1008    way_=-1;      // Swap direction
1009  }
1010  return 0.0;
1011}
1012// Print what would happen 
1013void
1014CbcIntegerBranchingObject::print()
1015{
1016  int iColumn = originalCbcObject_->columnNumber();
1017  assert (variable_==iColumn);
1018  if (way_<0) {
1019  { double olb,oub ;
1020    olb = model_->solver()->getColLower()[iColumn] ;
1021    oub = model_->solver()->getColUpper()[iColumn] ;
1022    printf("CbcInteger would branch down on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1023           iColumn,variable_,olb,oub,down_[0],down_[1]) ; }
1024  } else {
1025  { double olb,oub ;
1026    olb = model_->solver()->getColLower()[iColumn] ;
1027    oub = model_->solver()->getColUpper()[iColumn] ;
1028    printf("CbcInteger would branch up on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1029           iColumn,variable_,olb,oub,up_[0],up_[1]) ; }
1030  }
1031}
1032
1033
1034/** Default Constructor
1035
1036  Equivalent to an unspecified binary variable.
1037*/
1038CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ()
1039  : CbcSimpleInteger(),
1040    downPseudoCost_(1.0e-5),
1041    upPseudoCost_(1.0e-5),
1042    upDownSeparator_(-1.0),
1043    method_(0)
1044{
1045}
1046
1047/** Useful constructor
1048
1049  Loads actual upper & lower bounds for the specified variable.
1050*/
1051CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1052                                    int iColumn, double breakEven)
1053  : CbcSimpleInteger(model,iColumn,breakEven)
1054{
1055  const double * cost = model->getObjCoefficients();
1056  double costValue = CoinMax(1.0e-5,fabs(cost[iColumn]));
1057  // treat as if will cost what it says up
1058  upPseudoCost_=costValue;
1059  // and balance at breakeven
1060  downPseudoCost_=((1.0-breakEven_)*upPseudoCost_)/breakEven_;
1061  upDownSeparator_ = -1.0;
1062  method_=0;
1063}
1064
1065/** Useful constructor
1066
1067  Loads actual upper & lower bounds for the specified variable.
1068*/
1069CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1070                                    int iColumn, double downPseudoCost,
1071                                                        double upPseudoCost)
1072  : CbcSimpleInteger(model,iColumn)
1073{
1074  downPseudoCost_ = CoinMax(1.0e-10,downPseudoCost);
1075  upPseudoCost_ = CoinMax(1.0e-10,upPseudoCost);
1076  breakEven_ = upPseudoCost_/(upPseudoCost_+downPseudoCost_);
1077  upDownSeparator_ = -1.0;
1078  method_=0;
1079}
1080// Useful constructor - passed and model index and pseudo costs
1081CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model, int dummy,int iColumn, 
1082                                                        double downPseudoCost, double upPseudoCost)
1083{
1084  CbcSimpleIntegerPseudoCost(model,iColumn,downPseudoCost,upPseudoCost);
1085}
1086
1087// Copy constructor
1088CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ( const CbcSimpleIntegerPseudoCost & rhs)
1089  :CbcSimpleInteger(rhs),
1090   downPseudoCost_(rhs.downPseudoCost_),
1091   upPseudoCost_(rhs.upPseudoCost_),
1092   upDownSeparator_(rhs.upDownSeparator_),
1093   method_(rhs.method_)
1094
1095{
1096}
1097
1098// Clone
1099CbcObject *
1100CbcSimpleIntegerPseudoCost::clone() const
1101{
1102  return new CbcSimpleIntegerPseudoCost(*this);
1103}
1104
1105// Assignment operator
1106CbcSimpleIntegerPseudoCost & 
1107CbcSimpleIntegerPseudoCost::operator=( const CbcSimpleIntegerPseudoCost& rhs)
1108{
1109  if (this!=&rhs) {
1110    CbcSimpleInteger::operator=(rhs);
1111    downPseudoCost_=rhs.downPseudoCost_;
1112    upPseudoCost_=rhs.upPseudoCost_;
1113    upDownSeparator_=rhs.upDownSeparator_;
1114    method_=rhs.method_;
1115  }
1116  return *this;
1117}
1118
1119// Destructor
1120CbcSimpleIntegerPseudoCost::~CbcSimpleIntegerPseudoCost ()
1121{
1122}
1123// Creates a branching object
1124CbcBranchingObject * 
1125CbcSimpleIntegerPseudoCost::createBranch(int way) 
1126{
1127  OsiSolverInterface * solver = model_->solver();
1128  const double * solution = model_->testSolution();
1129  const double * lower = solver->getColLower();
1130  const double * upper = solver->getColUpper();
1131  double value = solution[columnNumber_];
1132  value = CoinMax(value, lower[columnNumber_]);
1133  value = CoinMin(value, upper[columnNumber_]);
1134#ifndef NDEBUG
1135  double nearest = floor(value+0.5);
1136  double integerTolerance = 
1137    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1138  assert (upper[columnNumber_]>lower[columnNumber_]);
1139#endif
1140  if (!model_->hotstartSolution()) {
1141    assert (fabs(value-nearest)>integerTolerance);
1142  } else {
1143    const double * hotstartSolution = model_->hotstartSolution();
1144    double targetValue = hotstartSolution[columnNumber_];
1145    if (way>0)
1146      value = targetValue-0.1;
1147    else
1148      value = targetValue+0.1;
1149  }
1150  CbcIntegerPseudoCostBranchingObject * newObject = 
1151    new CbcIntegerPseudoCostBranchingObject(model_,columnNumber_,way,
1152                                            value);
1153  double up =  upPseudoCost_*(ceil(value)-value);
1154  double down =  downPseudoCost_*(value-floor(value));
1155  double changeInGuessed=up-down;
1156  if (way>0)
1157    changeInGuessed = - changeInGuessed;
1158  changeInGuessed=CoinMax(0.0,changeInGuessed);
1159  //if (way>0)
1160  //changeInGuessed += 1.0e8; // bias to stay up
1161  newObject->setChangeInGuessed(changeInGuessed);
1162  newObject->setOriginalObject(this);
1163  return newObject;
1164}
1165// Infeasibility - large is 0.5
1166double 
1167CbcSimpleIntegerPseudoCost::infeasibility(int & preferredWay) const
1168{
1169  OsiSolverInterface * solver = model_->solver();
1170  const double * solution = model_->testSolution();
1171  const double * lower = solver->getColLower();
1172  const double * upper = solver->getColUpper();
1173  if (upper[columnNumber_]==lower[columnNumber_]) {
1174    // fixed
1175    preferredWay=1;
1176    return 0.0;
1177  }
1178  double value = solution[columnNumber_];
1179  value = CoinMax(value, lower[columnNumber_]);
1180  value = CoinMin(value, upper[columnNumber_]);
1181  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
1182    solution[columnNumber_],upper[columnNumber_]);*/
1183  double nearest = floor(value+0.5);
1184  double integerTolerance = 
1185    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1186  double below = floor(value+integerTolerance);
1187  double above = below+1.0;
1188  if (above>upper[columnNumber_]) {
1189    above=below;
1190    below = above -1;
1191  }
1192  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1193  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1194  if (downCost>=upCost)
1195    preferredWay=1;
1196  else
1197    preferredWay=-1;
1198  // See if up down choice set
1199  if (upDownSeparator_>0.0) {
1200    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
1201  }
1202  if (preferredWay_)
1203    preferredWay=preferredWay_;
1204  if (fabs(value-nearest)<=integerTolerance) {
1205    return 0.0;
1206  } else {
1207    // can't get at model so 1,2 don't make sense
1208    assert(method_<1||method_>2);
1209    if (!method_)
1210      return CoinMin(downCost,upCost);
1211    else
1212      return CoinMax(downCost,upCost);
1213  }
1214}
1215
1216// Return "up" estimate
1217double 
1218CbcSimpleIntegerPseudoCost::upEstimate() const
1219{
1220  OsiSolverInterface * solver = model_->solver();
1221  const double * solution = model_->testSolution();
1222  const double * lower = solver->getColLower();
1223  const double * upper = solver->getColUpper();
1224  double value = solution[columnNumber_];
1225  value = CoinMax(value, lower[columnNumber_]);
1226  value = CoinMin(value, upper[columnNumber_]);
1227  if (upper[columnNumber_]==lower[columnNumber_]) {
1228    // fixed
1229    return 0.0;
1230  }
1231  double integerTolerance = 
1232    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1233  double below = floor(value+integerTolerance);
1234  double above = below+1.0;
1235  if (above>upper[columnNumber_]) {
1236    above=below;
1237    below = above -1;
1238  }
1239  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1240  return upCost;
1241}
1242// Return "down" estimate
1243double 
1244CbcSimpleIntegerPseudoCost::downEstimate() const
1245{
1246  OsiSolverInterface * solver = model_->solver();
1247  const double * solution = model_->testSolution();
1248  const double * lower = solver->getColLower();
1249  const double * upper = solver->getColUpper();
1250  double value = solution[columnNumber_];
1251  value = CoinMax(value, lower[columnNumber_]);
1252  value = CoinMin(value, upper[columnNumber_]);
1253  if (upper[columnNumber_]==lower[columnNumber_]) {
1254    // fixed
1255    return 0.0;
1256  }
1257  double integerTolerance = 
1258    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1259  double below = floor(value+integerTolerance);
1260  double above = below+1.0;
1261  if (above>upper[columnNumber_]) {
1262    above=below;
1263    below = above -1;
1264  }
1265  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1266  return downCost;
1267}
1268
1269// Default Constructor
1270CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1271  :CbcIntegerBranchingObject()
1272{
1273  changeInGuessed_=1.0e-5;
1274}
1275
1276// Useful constructor
1277CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1278                                                      int variable, int way , double value)
1279  :CbcIntegerBranchingObject(model,variable,way,value)
1280{
1281}
1282// Useful constructor for fixing
1283CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1284                                                      int variable, int way,
1285                                                      double lowerValue, 
1286                                                      double upperValue)
1287  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
1288{
1289  changeInGuessed_=1.0e100;
1290}
1291 
1292
1293// Copy constructor
1294CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject ( 
1295                                 const CbcIntegerPseudoCostBranchingObject & rhs)
1296  :CbcIntegerBranchingObject(rhs)
1297{
1298  changeInGuessed_ = rhs.changeInGuessed_;
1299}
1300
1301// Assignment operator
1302CbcIntegerPseudoCostBranchingObject & 
1303CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject& rhs)
1304{
1305  if (this != &rhs) {
1306    CbcIntegerBranchingObject::operator=(rhs);
1307    changeInGuessed_ = rhs.changeInGuessed_;
1308  }
1309  return *this;
1310}
1311CbcBranchingObject * 
1312CbcIntegerPseudoCostBranchingObject::clone() const
1313{ 
1314  return (new CbcIntegerPseudoCostBranchingObject(*this));
1315}
1316
1317
1318// Destructor
1319CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1320{
1321}
1322
1323/*
1324  Perform a branch by adjusting the bounds of the specified variable. Note
1325  that each arm of the branch advances the object to the next arm by
1326  advancing the value of way_.
1327
1328  Providing new values for the variable's lower and upper bounds for each
1329  branching direction gives a little bit of additional flexibility and will
1330  be easily extensible to multi-way branching.
1331  Returns change in guessed objective on next branch
1332*/
1333double
1334CbcIntegerPseudoCostBranchingObject::branch()
1335{
1336  CbcIntegerBranchingObject::branch();
1337  return changeInGuessed_;
1338}
1339
1340
1341// Default Constructor
1342CbcCliqueBranchingObject::CbcCliqueBranchingObject()
1343  :CbcBranchingObject()
1344{
1345  clique_ = NULL;
1346  downMask_[0]=0;
1347  downMask_[1]=0;
1348  upMask_[0]=0;
1349  upMask_[1]=0;
1350}
1351
1352// Useful constructor
1353CbcCliqueBranchingObject::CbcCliqueBranchingObject (CbcModel * model,
1354                                                    const CbcClique * clique,
1355                                                    int way ,
1356                                                    int numberOnDownSide, const int * down,
1357                                                    int numberOnUpSide, const int * up)
1358  :CbcBranchingObject(model,clique->id(),way,0.5)
1359{
1360  clique_ = clique;
1361  downMask_[0]=0;
1362  downMask_[1]=0;
1363  upMask_[0]=0;
1364  upMask_[1]=0;
1365  int i;
1366  for (i=0;i<numberOnDownSide;i++) {
1367    int sequence = down[i];
1368    int iWord = sequence>>5;
1369    int iBit = sequence - 32*iWord;
1370    unsigned int k = 1<<iBit;
1371    downMask_[iWord] |= k;
1372  }
1373  for (i=0;i<numberOnUpSide;i++) {
1374    int sequence = up[i];
1375    int iWord = sequence>>5;
1376    int iBit = sequence - 32*iWord;
1377    unsigned int k = 1<<iBit;
1378    upMask_[iWord] |= k;
1379  }
1380}
1381
1382// Copy constructor
1383CbcCliqueBranchingObject::CbcCliqueBranchingObject ( const CbcCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1384{
1385  clique_=rhs.clique_;
1386  downMask_[0]=rhs.downMask_[0];
1387  downMask_[1]=rhs.downMask_[1];
1388  upMask_[0]=rhs.upMask_[0];
1389  upMask_[1]=rhs.upMask_[1];
1390}
1391
1392// Assignment operator
1393CbcCliqueBranchingObject & 
1394CbcCliqueBranchingObject::operator=( const CbcCliqueBranchingObject& rhs)
1395{
1396  if (this != &rhs) {
1397    CbcBranchingObject::operator=(rhs);
1398    clique_=rhs.clique_;
1399    downMask_[0]=rhs.downMask_[0];
1400    downMask_[1]=rhs.downMask_[1];
1401    upMask_[0]=rhs.upMask_[0];
1402    upMask_[1]=rhs.upMask_[1];
1403  }
1404  return *this;
1405}
1406CbcBranchingObject * 
1407CbcCliqueBranchingObject::clone() const
1408{ 
1409  return (new CbcCliqueBranchingObject(*this));
1410}
1411
1412
1413// Destructor
1414CbcCliqueBranchingObject::~CbcCliqueBranchingObject ()
1415{
1416}
1417double
1418CbcCliqueBranchingObject::branch()
1419{
1420  decrementNumberBranchesLeft();
1421  int iWord;
1422  int numberMembers = clique_->numberMembers();
1423  const int * which = clique_->members();
1424  const int * integerVariables = model_->integerVariable();
1425  int numberWords=(numberMembers+31)>>5;
1426  // *** for way - up means fix all those in down section
1427  if (way_<0) {
1428#ifdef FULL_PRINT
1429    printf("Down Fix ");
1430#endif
1431    for (iWord=0;iWord<numberWords;iWord++) {
1432      int i;
1433      for (i=0;i<32;i++) {
1434        unsigned int k = 1<<i;
1435        if ((upMask_[iWord]&k)!=0) {
1436          int iColumn = which[i+32*iWord];
1437#ifdef FULL_PRINT
1438          printf("%d ",i+32*iWord);
1439#endif
1440          // fix weak way
1441          if (clique_->type(i+32*iWord))
1442            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1443          else
1444            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1445        }
1446      }
1447    }
1448    way_=1;       // Swap direction
1449  } else {
1450#ifdef FULL_PRINT
1451    printf("Up Fix ");
1452#endif
1453    for (iWord=0;iWord<numberWords;iWord++) {
1454      int i;
1455      for (i=0;i<32;i++) {
1456        unsigned int k = 1<<i;
1457        if ((downMask_[iWord]&k)!=0) {
1458          int iColumn = which[i+32*iWord];
1459#ifdef FULL_PRINT
1460          printf("%d ",i+32*iWord);
1461#endif
1462          // fix weak way
1463          if (clique_->type(i+32*iWord))
1464            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1465          else
1466            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1467        }
1468      }
1469    }
1470    way_=-1;      // Swap direction
1471  }
1472#ifdef FULL_PRINT
1473  printf("\n");
1474#endif
1475  return 0.0;
1476}
1477// Print what would happen 
1478void
1479CbcCliqueBranchingObject::print()
1480{
1481  int iWord;
1482  int numberMembers = clique_->numberMembers();
1483  const int * which = clique_->members();
1484  const int * integerVariables = model_->integerVariable();
1485  int numberWords=(numberMembers+31)>>5;
1486  // *** for way - up means fix all those in down section
1487  if (way_<0) {
1488    printf("Clique - Down Fix ");
1489    for (iWord=0;iWord<numberWords;iWord++) {
1490      int i;
1491      for (i=0;i<32;i++) {
1492        unsigned int k = 1<<i;
1493        if ((upMask_[iWord]&k)!=0) {
1494          int iColumn = which[i+32*iWord];
1495          printf("%d ",integerVariables[iColumn]);
1496        }
1497      }
1498    }
1499  } else {
1500    printf("Clique - Up Fix ");
1501    for (iWord=0;iWord<numberWords;iWord++) {
1502      int i;
1503      for (i=0;i<32;i++) {
1504        unsigned int k = 1<<i;
1505        if ((downMask_[iWord]&k)!=0) {
1506          int iColumn = which[i+32*iWord];
1507          printf("%d ",integerVariables[iColumn]);
1508        }
1509      }
1510    }
1511  }
1512  printf("\n");
1513}
1514 
1515// Default Constructor
1516CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject()
1517  :CbcBranchingObject()
1518{
1519  clique_=NULL;
1520  downMask_=NULL;
1521  upMask_=NULL;
1522}
1523
1524// Useful constructor
1525CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject (CbcModel * model,
1526                                                            const CbcClique * clique, 
1527                                                            int way ,
1528                                                    int numberOnDownSide, const int * down,
1529                                                    int numberOnUpSide, const int * up)
1530  :CbcBranchingObject(model,clique->id(),way,0.5)
1531{
1532  clique_ = clique;
1533  int numberMembers = clique_->numberMembers();
1534  int numberWords=(numberMembers+31)>>5;
1535  downMask_ = new unsigned int [numberWords];
1536  upMask_ = new unsigned int [numberWords];
1537  memset(downMask_,0,numberWords*sizeof(unsigned int));
1538  memset(upMask_,0,numberWords*sizeof(unsigned int));
1539  int i;
1540  for (i=0;i<numberOnDownSide;i++) {
1541    int sequence = down[i];
1542    int iWord = sequence>>5;
1543    int iBit = sequence - 32*iWord;
1544    unsigned int k = 1<<iBit;
1545    downMask_[iWord] |= k;
1546  }
1547  for (i=0;i<numberOnUpSide;i++) {
1548    int sequence = up[i];
1549    int iWord = sequence>>5;
1550    int iBit = sequence - 32*iWord;
1551    unsigned int k = 1<<iBit;
1552    upMask_[iWord] |= k;
1553  }
1554}
1555
1556// Copy constructor
1557CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject ( const CbcLongCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1558{
1559  clique_=rhs.clique_;
1560  if (rhs.downMask_) {
1561    int numberMembers = clique_->numberMembers();
1562    int numberWords=(numberMembers+31)>>5;
1563    downMask_ = new unsigned int [numberWords];
1564    memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
1565    upMask_ = new unsigned int [numberWords];
1566    memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
1567  } else {
1568    downMask_=NULL;
1569    upMask_=NULL;
1570  }   
1571}
1572
1573// Assignment operator
1574CbcLongCliqueBranchingObject & 
1575CbcLongCliqueBranchingObject::operator=( const CbcLongCliqueBranchingObject& rhs)
1576{
1577  if (this != &rhs) {
1578    CbcBranchingObject::operator=(rhs);
1579    clique_=rhs.clique_;
1580    delete [] downMask_;
1581    delete [] upMask_;
1582    if (rhs.downMask_) {
1583      int numberMembers = clique_->numberMembers();
1584      int numberWords=(numberMembers+31)>>5;
1585      downMask_ = new unsigned int [numberWords];
1586      memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
1587      upMask_ = new unsigned int [numberWords];
1588      memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
1589    } else {
1590      downMask_=NULL;
1591      upMask_=NULL;
1592    }   
1593  }
1594  return *this;
1595}
1596CbcBranchingObject * 
1597CbcLongCliqueBranchingObject::clone() const
1598{ 
1599  return (new CbcLongCliqueBranchingObject(*this));
1600}
1601
1602
1603// Destructor
1604CbcLongCliqueBranchingObject::~CbcLongCliqueBranchingObject ()
1605{
1606  delete [] downMask_;
1607  delete [] upMask_;
1608}
1609double
1610CbcLongCliqueBranchingObject::branch()
1611{
1612  decrementNumberBranchesLeft();
1613  int iWord;
1614  int numberMembers = clique_->numberMembers();
1615  const int * which = clique_->members();
1616  const int * integerVariables = model_->integerVariable();
1617  int numberWords=(numberMembers+31)>>5;
1618  // *** for way - up means fix all those in down section
1619  if (way_<0) {
1620#ifdef FULL_PRINT
1621    printf("Down Fix ");
1622#endif
1623    for (iWord=0;iWord<numberWords;iWord++) {
1624      int i;
1625      for (i=0;i<32;i++) {
1626        unsigned int k = 1<<i;
1627        if ((upMask_[iWord]&k)!=0) {
1628          int iColumn = which[i+32*iWord];
1629#ifdef FULL_PRINT
1630          printf("%d ",i+32*iWord);
1631#endif
1632          // fix weak way
1633          if (clique_->type(i+32*iWord))
1634            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1635          else
1636            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1637        }
1638      }
1639    }
1640    way_=1;       // Swap direction
1641  } else {
1642#ifdef FULL_PRINT
1643    printf("Up Fix ");
1644#endif
1645    for (iWord=0;iWord<numberWords;iWord++) {
1646      int i;
1647      for (i=0;i<32;i++) {
1648        unsigned int k = 1<<i;
1649        if ((downMask_[iWord]&k)!=0) {
1650          int iColumn = which[i+32*iWord];
1651#ifdef FULL_PRINT
1652          printf("%d ",i+32*iWord);
1653#endif
1654          // fix weak way
1655          if (clique_->type(i+32*iWord))
1656            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1657          else
1658            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1659        }
1660      }
1661    }
1662    way_=-1;      // Swap direction
1663  }
1664#ifdef FULL_PRINT
1665  printf("\n");
1666#endif
1667  return 0.0;
1668}
1669void
1670CbcLongCliqueBranchingObject::print()
1671{
1672  int iWord;
1673  int numberMembers = clique_->numberMembers();
1674  const int * which = clique_->members();
1675  const int * integerVariables = model_->integerVariable();
1676  int numberWords=(numberMembers+31)>>5;
1677  // *** for way - up means fix all those in down section
1678  if (way_<0) {
1679    printf("Clique - Down Fix ");
1680    for (iWord=0;iWord<numberWords;iWord++) {
1681      int i;
1682      for (i=0;i<32;i++) {
1683        unsigned int k = 1<<i;
1684        if ((upMask_[iWord]&k)!=0) {
1685          int iColumn = which[i+32*iWord];
1686          printf("%d ",integerVariables[iColumn]);
1687        }
1688      }
1689    }
1690  } else {
1691    printf("Clique - Up Fix ");
1692    for (iWord=0;iWord<numberWords;iWord++) {
1693      int i;
1694      for (i=0;i<32;i++) {
1695        unsigned int k = 1<<i;
1696        if ((downMask_[iWord]&k)!=0) {
1697          int iColumn = which[i+32*iWord];
1698          printf("%d ",integerVariables[iColumn]);
1699        }
1700      }
1701    }
1702  }
1703  printf("\n");
1704}
1705// Default Constructor
1706CbcSOSBranchingObject::CbcSOSBranchingObject()
1707  :CbcBranchingObject()
1708{
1709  set_ = NULL;
1710  separator_=0.0;
1711}
1712
1713// Useful constructor
1714CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
1715                                              const CbcSOS * set,
1716                                              int way ,
1717                                              double separator)
1718  :CbcBranchingObject(model,set->id(),way,0.5)
1719{
1720  set_ = set;
1721  separator_ = separator;
1722}
1723
1724// Copy constructor
1725CbcSOSBranchingObject::CbcSOSBranchingObject ( const CbcSOSBranchingObject & rhs) :CbcBranchingObject(rhs)
1726{
1727  set_=rhs.set_;
1728  separator_ = rhs.separator_;
1729}
1730
1731// Assignment operator
1732CbcSOSBranchingObject & 
1733CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject& rhs)
1734{
1735  if (this != &rhs) {
1736    CbcBranchingObject::operator=(rhs);
1737    set_=rhs.set_;
1738    separator_ = rhs.separator_;
1739  }
1740  return *this;
1741}
1742CbcBranchingObject * 
1743CbcSOSBranchingObject::clone() const
1744{ 
1745  return (new CbcSOSBranchingObject(*this));
1746}
1747
1748
1749// Destructor
1750CbcSOSBranchingObject::~CbcSOSBranchingObject ()
1751{
1752}
1753double
1754CbcSOSBranchingObject::branch()
1755{
1756  decrementNumberBranchesLeft();
1757  int numberMembers = set_->numberMembers();
1758  const int * which = set_->members();
1759  const double * weights = set_->weights();
1760  OsiSolverInterface * solver = model_->solver();
1761  //const double * lower = solver->getColLower();
1762  //const double * upper = solver->getColUpper();
1763  // *** for way - up means fix all those in down section
1764  if (way_<0) {
1765    int i;
1766    for ( i=0;i<numberMembers;i++) {
1767      if (weights[i] > separator_)
1768        break;
1769    }
1770    assert (i<numberMembers);
1771    for (;i<numberMembers;i++) 
1772      solver->setColUpper(which[i],0.0);
1773    way_=1;       // Swap direction
1774  } else {
1775    int i;
1776    for ( i=0;i<numberMembers;i++) {
1777      if (weights[i] >= separator_)
1778        break;
1779      else
1780        solver->setColUpper(which[i],0.0);
1781    }
1782    assert (i<numberMembers);
1783    way_=-1;      // Swap direction
1784  }
1785  return 0.0;
1786}
1787// Print what would happen 
1788void
1789CbcSOSBranchingObject::print()
1790{
1791  int numberMembers = set_->numberMembers();
1792  const int * which = set_->members();
1793  const double * weights = set_->weights();
1794  OsiSolverInterface * solver = model_->solver();
1795  //const double * lower = solver->getColLower();
1796  const double * upper = solver->getColUpper();
1797  int first=numberMembers;
1798  int last=-1;
1799  int numberFixed=0;
1800  int numberOther=0;
1801  int i;
1802  for ( i=0;i<numberMembers;i++) {
1803    double bound = upper[which[i]];
1804    if (bound) {
1805      first = CoinMin(first,i);
1806      last = CoinMax(last,i);
1807    }
1808  }
1809  // *** for way - up means fix all those in down section
1810  if (way_<0) {
1811    printf("SOS Down");
1812    for ( i=0;i<numberMembers;i++) {
1813      double bound = upper[which[i]];
1814      if (weights[i] > separator_)
1815        break;
1816      else if (bound)
1817        numberOther++;
1818    }
1819    assert (i<numberMembers);
1820    for (;i<numberMembers;i++) {
1821      double bound = upper[which[i]];
1822      if (bound)
1823        numberFixed++;
1824    }
1825  } else {
1826    printf("SOS Up");
1827    for ( i=0;i<numberMembers;i++) {
1828      double bound = upper[which[i]];
1829      if (weights[i] >= separator_)
1830        break;
1831      else if (bound)
1832        numberFixed++;
1833    }
1834    assert (i<numberMembers);
1835    for (;i<numberMembers;i++) {
1836      double bound = upper[which[i]];
1837      if (bound)
1838        numberOther++;
1839    }
1840  }
1841  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
1842         separator_,which[first],weights[first],which[last],weights[last],numberFixed,numberOther);
1843}
1844 
1845// Default Constructor
1846CbcBranchDefaultDecision::CbcBranchDefaultDecision()
1847  :CbcBranchDecision()
1848{
1849  bestCriterion_ = 0.0;
1850  bestChangeUp_ = 0.0;
1851  bestNumberUp_ = 0;
1852  bestChangeDown_ = 0.0;
1853  bestObject_ = NULL;
1854  bestNumberDown_ = 0;
1855}
1856
1857// Copy constructor
1858CbcBranchDefaultDecision::CbcBranchDefaultDecision (
1859                                    const CbcBranchDefaultDecision & rhs)
1860  :CbcBranchDecision(rhs)
1861{
1862  bestCriterion_ = rhs.bestCriterion_;
1863  bestChangeUp_ = rhs.bestChangeUp_;
1864  bestNumberUp_ = rhs.bestNumberUp_;
1865  bestChangeDown_ = rhs.bestChangeDown_;
1866  bestNumberDown_ = rhs.bestNumberDown_;
1867  bestObject_ = rhs.bestObject_;
1868  model_ = rhs.model_;
1869}
1870
1871CbcBranchDefaultDecision::~CbcBranchDefaultDecision()
1872{
1873}
1874
1875// Clone
1876CbcBranchDecision * 
1877CbcBranchDefaultDecision::clone() const
1878{
1879  return new CbcBranchDefaultDecision(*this);
1880}
1881
1882// Initialize i.e. before start of choosing at a node
1883void 
1884CbcBranchDefaultDecision::initialize(CbcModel * model)
1885{
1886  bestCriterion_ = 0.0;
1887  bestChangeUp_ = 0.0;
1888  bestNumberUp_ = 0;
1889  bestChangeDown_ = 0.0;
1890  bestNumberDown_ = 0;
1891  bestObject_ = NULL;
1892  model_ = model;
1893}
1894
1895
1896/*
1897  Simple default decision algorithm. Compare based on infeasibility (numInfUp,
1898  numInfDn) until a solution is found by search, then switch to change in
1899  objective (changeUp, changeDn). Note that bestSoFar is remembered in
1900  bestObject_, so the parameter bestSoFar is unused.
1901*/
1902
1903int
1904CbcBranchDefaultDecision::betterBranch(CbcBranchingObject * thisOne,
1905                            CbcBranchingObject * bestSoFar,
1906                            double changeUp, int numInfUp,
1907                            double changeDn, int numInfDn)
1908{
1909  bool beforeSolution = cbcModel()->getSolutionCount()==
1910    cbcModel()->getNumberHeuristicSolutions();;
1911  int betterWay=0;
1912  if (beforeSolution) {
1913    if (!bestObject_) {
1914      bestNumberUp_=INT_MAX;
1915      bestNumberDown_=INT_MAX;
1916    }
1917    // before solution - choose smallest number
1918    // could add in depth as well
1919    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1920    if (numInfUp<numInfDn) {
1921      if (numInfUp<bestNumber) {
1922        betterWay = 1;
1923      } else if (numInfUp==bestNumber) {
1924        if (changeUp<bestCriterion_)
1925          betterWay=1;
1926      }
1927    } else if (numInfUp>numInfDn) {
1928      if (numInfDn<bestNumber) {
1929        betterWay = -1;
1930      } else if (numInfDn==bestNumber) {
1931        if (changeDn<bestCriterion_)
1932          betterWay=-1;
1933      }
1934    } else {
1935      // up and down have same number
1936      bool better=false;
1937      if (numInfUp<bestNumber) {
1938        better=true;
1939      } else if (numInfUp==bestNumber) {
1940        if (min(changeUp,changeDn)<bestCriterion_)
1941          better=true;;
1942      }
1943      if (better) {
1944        // see which way
1945        if (changeUp<=changeDn)
1946          betterWay=1;
1947        else
1948          betterWay=-1;
1949      }
1950    }
1951  } else {
1952    if (!bestObject_) {
1953      bestCriterion_=-1.0;
1954    }
1955    // got a solution
1956    if (changeUp<=changeDn) {
1957      if (changeUp>bestCriterion_)
1958        betterWay=1;
1959    } else {
1960      if (changeDn>bestCriterion_)
1961        betterWay=-1;
1962    }
1963  }
1964  if (betterWay) {
1965    bestCriterion_ = CoinMin(changeUp,changeDn);
1966    bestChangeUp_ = changeUp;
1967    bestNumberUp_ = numInfUp;
1968    bestChangeDown_ = changeDn;
1969    bestNumberDown_ = numInfDn;
1970    bestObject_=thisOne;
1971    // See if user is overriding way
1972    if (thisOne->object()&&thisOne->object()->preferredWay())
1973      betterWay = thisOne->object()->preferredWay();
1974  }
1975  return betterWay;
1976}
1977/* Sets or gets best criterion so far */
1978void 
1979CbcBranchDefaultDecision::setBestCriterion(double value)
1980{ 
1981  bestCriterion_ = value;
1982}
1983double 
1984CbcBranchDefaultDecision::getBestCriterion() const
1985{ 
1986  return bestCriterion_;
1987}
1988
1989/* Compare N branching objects. Return index of best
1990   and sets way of branching in chosen object.
1991   
1992   This routine is used only after strong branching.
1993*/
1994
1995int
1996CbcBranchDefaultDecision::bestBranch (CbcBranchingObject ** objects, int numberObjects,
1997                                   int numberUnsatisfied,
1998                                   double * changeUp, int * numberInfeasibilitiesUp,
1999                                   double * changeDown, int * numberInfeasibilitiesDown,
2000                                   double objectiveValue) 
2001{
2002
2003  int bestWay=0;
2004  int whichObject = -1;
2005  if (numberObjects) {
2006    CbcModel * model = cbcModel();
2007    // at continuous
2008    //double continuousObjective = model->getContinuousObjective();
2009    //int continuousInfeasibilities = model->getContinuousInfeasibilities();
2010   
2011    // average cost to get rid of infeasibility
2012    //double averageCostPerInfeasibility =
2013    //(objectiveValue-continuousObjective)/
2014    //(double) (abs(continuousInfeasibilities-numberUnsatisfied)+1);
2015    /* beforeSolution is :
2016       0 - before any solution
2017       n - n heuristic solutions but no branched one
2018       -1 - branched solution found
2019    */
2020    int numberSolutions = model->getSolutionCount();
2021    double cutoff = model->getCutoff();
2022    int method=0;
2023    int i;
2024    if (numberSolutions) {
2025      int numberHeuristic = model->getNumberHeuristicSolutions();
2026      if (numberHeuristic<numberSolutions) {
2027        method = 1;
2028      } else {
2029        method = 2;
2030        // look further
2031        for ( i = 0 ; i < numberObjects ; i++) {
2032          int numberNext = numberInfeasibilitiesUp[i];
2033         
2034          if (numberNext<numberUnsatisfied) {
2035            int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
2036            double perUnsatisfied = changeUp[i]/(double) numberUp;
2037            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2038            if (estimatedObjective<cutoff) 
2039              method=3;
2040          }
2041          numberNext = numberInfeasibilitiesDown[i];
2042          if (numberNext<numberUnsatisfied) {
2043            int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
2044            double perUnsatisfied = changeDown[i]/(double) numberDown;
2045            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2046            if (estimatedObjective<cutoff) 
2047              method=3;
2048          }
2049        }
2050      }
2051      method=2;
2052    } else {
2053      method = 0;
2054    }
2055    // Uncomment next to force method 4
2056    //method=4;
2057    /* Methods :
2058       0 - fewest infeasibilities
2059       1 - largest min change in objective
2060       2 - as 1 but use sum of changes if min close
2061       3 - predicted best solution
2062       4 - take cheapest up branch if infeasibilities same
2063    */
2064    int bestNumber=INT_MAX;
2065    double bestCriterion=-1.0e50;
2066    double alternativeCriterion = -1.0;
2067    double bestEstimate = 1.0e100;
2068    switch (method) {
2069    case 0:
2070      // could add in depth as well
2071      for ( i = 0 ; i < numberObjects ; i++) {
2072        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
2073        if (thisNumber<=bestNumber) {
2074          int betterWay=0;
2075          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
2076            if (numberInfeasibilitiesUp[i]<bestNumber) {
2077              betterWay = 1;
2078            } else {
2079              if (changeUp[i]<bestCriterion)
2080                betterWay=1;
2081            }
2082          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
2083            if (numberInfeasibilitiesDown[i]<bestNumber) {
2084              betterWay = -1;
2085            } else {
2086              if (changeDown[i]<bestCriterion)
2087                betterWay=-1;
2088            }
2089          } else {
2090            // up and down have same number
2091            bool better=false;
2092            if (numberInfeasibilitiesUp[i]<bestNumber) {
2093              better=true;
2094            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2095              if (min(changeUp[i],changeDown[i])<bestCriterion)
2096                better=true;;
2097            }
2098            if (better) {
2099              // see which way
2100              if (changeUp[i]<=changeDown[i])
2101                betterWay=1;
2102              else
2103                betterWay=-1;
2104            }
2105          }
2106          if (betterWay) {
2107            bestCriterion = min(changeUp[i],changeDown[i]);
2108            bestNumber = thisNumber;
2109            whichObject = i;
2110            bestWay = betterWay;
2111          }
2112        }
2113      }
2114      break;
2115    case 1:
2116      for ( i = 0 ; i < numberObjects ; i++) {
2117        int betterWay=0;
2118        if (changeUp[i]<=changeDown[i]) {
2119          if (changeUp[i]>bestCriterion)
2120            betterWay=1;
2121        } else {
2122          if (changeDown[i]>bestCriterion)
2123            betterWay=-1;
2124        }
2125        if (betterWay) {
2126          bestCriterion = min(changeUp[i],changeDown[i]);
2127          whichObject = i;
2128          bestWay = betterWay;
2129        }
2130      }
2131      break;
2132    case 2:
2133      for ( i = 0 ; i < numberObjects ; i++) {
2134        double change = min(changeUp[i],changeDown[i]);
2135        double sum = changeUp[i] + changeDown[i];
2136        bool take=false;
2137        if (change>1.1*bestCriterion) 
2138          take=true;
2139        else if (change>0.9*bestCriterion&&sum+change>bestCriterion+alternativeCriterion) 
2140          take=true;
2141        if (take) {
2142          if (changeUp[i]<=changeDown[i]) {
2143            if (changeUp[i]>bestCriterion)
2144              bestWay=1;
2145          } else {
2146            if (changeDown[i]>bestCriterion)
2147              bestWay=-1;
2148          }
2149          bestCriterion = change;
2150          alternativeCriterion = sum;
2151          whichObject = i;
2152        }
2153      }
2154      break;
2155    case 3:
2156      for ( i = 0 ; i < numberObjects ; i++) {
2157        int numberNext = numberInfeasibilitiesUp[i];
2158       
2159        if (numberNext<numberUnsatisfied) {
2160          int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
2161          double perUnsatisfied = changeUp[i]/(double) numberUp;
2162          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2163          if (estimatedObjective<bestEstimate) {
2164            bestEstimate = estimatedObjective;
2165            bestWay=1;
2166            whichObject=i;
2167          }
2168        }
2169        numberNext = numberInfeasibilitiesDown[i];
2170        if (numberNext<numberUnsatisfied) {
2171          int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
2172          double perUnsatisfied = changeDown[i]/(double) numberDown;
2173          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2174          if (estimatedObjective<bestEstimate) {
2175            bestEstimate = estimatedObjective;
2176            bestWay=-1;
2177            whichObject=i;
2178          }
2179        }
2180      }
2181      break;
2182    case 4:
2183      // if number infeas same then cheapest up
2184      // first get best number or when going down
2185      // now choose smallest change up amongst equal number infeas
2186      for ( i = 0 ; i < numberObjects ; i++) {
2187        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
2188        if (thisNumber<=bestNumber) {
2189          int betterWay=0;
2190          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
2191            if (numberInfeasibilitiesUp[i]<bestNumber) {
2192              betterWay = 1;
2193            } else {
2194              if (changeUp[i]<bestCriterion)
2195                betterWay=1;
2196            }
2197          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
2198            if (numberInfeasibilitiesDown[i]<bestNumber) {
2199              betterWay = -1;
2200            } else {
2201              if (changeDown[i]<bestCriterion)
2202                betterWay=-1;
2203            }
2204          } else {
2205            // up and down have same number
2206            bool better=false;
2207            if (numberInfeasibilitiesUp[i]<bestNumber) {
2208              better=true;
2209            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2210              if (min(changeUp[i],changeDown[i])<bestCriterion)
2211                better=true;;
2212            }
2213            if (better) {
2214              // see which way
2215              if (changeUp[i]<=changeDown[i])
2216                betterWay=1;
2217              else
2218                betterWay=-1;
2219            }
2220          }
2221          if (betterWay) {
2222            bestCriterion = min(changeUp[i],changeDown[i]);
2223            bestNumber = thisNumber;
2224            whichObject = i;
2225            bestWay = betterWay;
2226          }
2227        }
2228      }
2229      bestCriterion=1.0e50;
2230      for ( i = 0 ; i < numberObjects ; i++) {
2231        int thisNumber = numberInfeasibilitiesUp[i];
2232        if (thisNumber==bestNumber&&changeUp) {
2233          if (changeUp[i]<bestCriterion) {
2234            bestCriterion = changeUp[i];
2235            whichObject = i;
2236            bestWay = 1;
2237          }
2238        }
2239      }
2240      break;
2241    }
2242    // set way in best
2243    if (whichObject>=0) {
2244      CbcBranchingObject * bestObject = objects[whichObject];
2245      if (bestObject->object()&&bestObject->object()->preferredWay()) 
2246        bestWay = bestObject->object()->preferredWay();
2247      bestObject->way(bestWay);
2248    } else {
2249      printf("debug\n");
2250    }
2251  }
2252  return whichObject;
2253}
2254
2255// Default Constructor
2256CbcFollowOn::CbcFollowOn ()
2257  : CbcObject(),
2258    rhs_(NULL)
2259{
2260}
2261
2262// Useful constructor
2263CbcFollowOn::CbcFollowOn (CbcModel * model)
2264  : CbcObject(model)
2265{
2266  assert (model);
2267  OsiSolverInterface * solver = model_->solver();
2268  matrix_ = *solver->getMatrixByCol();
2269  matrix_.removeGaps();
2270  matrix_.setExtraGap(0.0);
2271  matrixByRow_ = *solver->getMatrixByRow();
2272  int numberRows = matrix_.getNumRows();
2273 
2274  rhs_ = new int[numberRows];
2275  int i;
2276  const double * rowLower = solver->getRowLower();
2277  const double * rowUpper = solver->getRowUpper();
2278  // Row copy
2279  const double * elementByRow = matrixByRow_.getElements();
2280  const int * column = matrixByRow_.getIndices();
2281  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2282  const int * rowLength = matrixByRow_.getVectorLengths();
2283  for (i=0;i<numberRows;i++) {
2284    rhs_[i]=0;
2285    double value = rowLower[i];
2286    if (value==rowUpper[i]) {
2287      if (floor(value)==value&&value>=1.0&&value<10.0) {
2288        // check elements
2289        bool good=true;
2290        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2291          int iColumn = column[j];
2292          if (!solver->isBinary(iColumn))
2293            good=false;
2294          double elValue = elementByRow[j];
2295          if (floor(elValue)!=elValue||value<1.0)
2296            good=false;
2297        }
2298        if (good)
2299          rhs_[i]=(int) value;
2300      }
2301    }
2302  }
2303}
2304
2305// Copy constructor
2306CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
2307  :CbcObject(rhs),
2308   matrix_(rhs.matrix_),
2309   matrixByRow_(rhs.matrixByRow_)
2310{
2311  int numberRows = matrix_.getNumRows();
2312  rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2313}
2314
2315// Clone
2316CbcObject *
2317CbcFollowOn::clone() const
2318{
2319  return new CbcFollowOn(*this);
2320}
2321
2322// Assignment operator
2323CbcFollowOn & 
2324CbcFollowOn::operator=( const CbcFollowOn& rhs)
2325{
2326  if (this!=&rhs) {
2327    CbcObject::operator=(rhs);
2328    delete [] rhs_;
2329    matrix_ = rhs.matrix_;
2330    matrixByRow_ = rhs.matrixByRow_;
2331    int numberRows = matrix_.getNumRows();
2332    rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2333  }
2334  return *this;
2335}
2336
2337// Destructor
2338CbcFollowOn::~CbcFollowOn ()
2339{
2340  delete [] rhs_;
2341}
2342// As some computation is needed in more than one place - returns row
2343int 
2344CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
2345{
2346  int whichRow=-1;
2347  otherRow=-1;
2348  int numberRows = matrix_.getNumRows();
2349 
2350  int i;
2351  // For sorting
2352  int * sort = new int [numberRows];
2353  int * isort = new int [numberRows];
2354  // Column copy
2355  //const double * element = matrix_.getElements();
2356  const int * row = matrix_.getIndices();
2357  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2358  const int * columnLength = matrix_.getVectorLengths();
2359  // Row copy
2360  const double * elementByRow = matrixByRow_.getElements();
2361  const int * column = matrixByRow_.getIndices();
2362  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2363  const int * rowLength = matrixByRow_.getVectorLengths();
2364  OsiSolverInterface * solver = model_->solver();
2365  const double * columnLower = solver->getColLower();
2366  const double * columnUpper = solver->getColUpper();
2367  const double * solution = solver->getColSolution();
2368  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2369  int nSort=0;
2370  for (i=0;i<numberRows;i++) {
2371    if (rhs_[i]) {
2372      // check elements
2373      double smallest=1.0e10;
2374      double largest=0.0;
2375      int rhsValue=rhs_[i];
2376      int number1=0;
2377      int numberUnsatisfied=0;
2378      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2379        int iColumn = column[j];
2380        double value = elementByRow[j];
2381        double solValue = solution[iColumn];
2382        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2383          smallest = CoinMin(smallest,value);
2384          largest = CoinMax(largest,value);
2385          if (value==1.0)
2386            number1++;
2387          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
2388            numberUnsatisfied++;
2389        } else {
2390          rhsValue -= (int)(value*floor(solValue+0.5));
2391        }
2392      }
2393      if (numberUnsatisfied>1) {
2394        if (smallest<largest) {
2395          // probably no good but check a few things
2396          assert (largest<=rhsValue);
2397          if (number1==1&&largest==rhsValue)
2398            printf("could fix\n");
2399        } else if (largest==rhsValue) {
2400          sort[nSort]=i;
2401          isort[nSort++]=-numberUnsatisfied;
2402        }
2403      }
2404    }
2405  }
2406  if (nSort>1) {
2407    CoinSort_2(isort,isort+nSort,sort);
2408    CoinZeroN(isort,numberRows);
2409    double * other = new double[numberRows];
2410    CoinZeroN(other,numberRows);
2411    int * which = new int[numberRows];
2412    //#define COUNT
2413#ifndef COUNT
2414    bool beforeSolution = model_->getSolutionCount()==0;
2415#endif
2416    for (int k=0;k<nSort-1;k++) {
2417      i=sort[k];
2418      int numberUnsatisfied = 0;
2419      int n=0;
2420      int j;
2421      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2422        int iColumn = column[j];
2423        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2424          double solValue = solution[iColumn]-columnLower[iColumn];
2425          if (solValue<1.0-integerTolerance&&solValue>integerTolerance) {
2426            numberUnsatisfied++;
2427            for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2428              int iRow = row[jj];
2429              if (rhs_[iRow]) {
2430                other[iRow]+=solValue;
2431                if (isort[iRow]) {
2432                  isort[iRow]++;
2433                } else {
2434                  isort[iRow]=1;
2435                  which[n++]=iRow;
2436                }
2437              }
2438            }
2439          }
2440        }
2441      }
2442      double total=0.0;
2443      // Take out row
2444      double sumThis=other[i];
2445      other[i]=0.0;
2446      assert (numberUnsatisfied==isort[i]);
2447      // find one nearest half if solution, one if before solution
2448      int iBest=-1;
2449      double dtarget=0.5*total;
2450#ifdef COUNT
2451      int target = (numberUnsatisfied+1)>>1;
2452      int best=numberUnsatisfied;
2453#else
2454      double best;
2455      if (beforeSolution)
2456        best=dtarget;
2457      else
2458        best=1.0e30;
2459#endif
2460      for (j=0;j<n;j++) {
2461        int iRow = which[j];
2462        double dvalue=other[iRow];
2463        other[iRow]=0.0;
2464#ifdef COUNT
2465        int value = isort[iRow];
2466#endif
2467        isort[iRow]=0;
2468        if (fabs(dvalue)<1.0e-8||fabs(sumThis-dvalue)<1.0e-8)
2469          continue;
2470        if (dvalue<integerTolerance||dvalue>1.0-integerTolerance)
2471          continue;
2472#ifdef COUNT
2473        if (abs(value-target)<best&&value!=numberUnsatisfied) {
2474          best=abs(value-target);
2475          iBest=iRow;
2476          if (dvalue<dtarget)
2477            preferredWay=1;
2478          else
2479            preferredWay=-1;
2480        }
2481#else
2482        if (beforeSolution) {
2483          if (fabs(dvalue-dtarget)>best) {
2484            best = fabs(dvalue-dtarget);
2485            iBest=iRow;
2486            if (dvalue<dtarget)
2487              preferredWay=1;
2488            else
2489              preferredWay=-1;
2490          }
2491        } else {
2492          if (fabs(dvalue-dtarget)<best) {
2493            best = fabs(dvalue-dtarget);
2494            iBest=iRow;
2495            if (dvalue<dtarget)
2496              preferredWay=1;
2497            else
2498              preferredWay=-1;
2499          }
2500        }
2501#endif
2502      }
2503      if (iBest>=0) {
2504        whichRow=i;
2505        otherRow=iBest;
2506        break;
2507      }
2508    }
2509    delete [] which;
2510    delete [] other;
2511  }
2512  delete [] sort;
2513  delete [] isort;
2514  return whichRow;
2515}
2516
2517// Infeasibility - large is 0.5
2518double 
2519CbcFollowOn::infeasibility(int & preferredWay) const
2520{
2521  int otherRow=0;
2522  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2523  if (whichRow<0)
2524    return 0.0;
2525  else
2526  return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
2527}
2528
2529// This looks at solution and sets bounds to contain solution
2530void 
2531CbcFollowOn::feasibleRegion()
2532{
2533}
2534
2535
2536// Creates a branching object
2537CbcBranchingObject * 
2538CbcFollowOn::createBranch(int way) 
2539{
2540  int otherRow=0;
2541  int preferredWay;
2542  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2543  assert(way==preferredWay);
2544  assert (whichRow>=0);
2545  int numberColumns = matrix_.getNumCols();
2546 
2547  // Column copy
2548  //const double * element = matrix_.getElements();
2549  const int * row = matrix_.getIndices();
2550  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2551  const int * columnLength = matrix_.getVectorLengths();
2552  // Row copy
2553  //const double * elementByRow = matrixByRow_.getElements();
2554  const int * column = matrixByRow_.getIndices();
2555  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2556  const int * rowLength = matrixByRow_.getVectorLengths();
2557  OsiSolverInterface * solver = model_->solver();
2558  const double * columnLower = solver->getColLower();
2559  const double * columnUpper = solver->getColUpper();
2560  //const double * solution = solver->getColSolution();
2561  int nUp=0;
2562  int nDown=0;
2563  int * upList = new int[numberColumns];
2564  int * downList = new int[numberColumns];
2565  int j;
2566  for (j=rowStart[whichRow];j<rowStart[whichRow]+rowLength[whichRow];j++) {
2567    int iColumn = column[j];
2568    if (columnLower[iColumn]!=columnUpper[iColumn]) {
2569      bool up=true;
2570      for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2571        int iRow = row[jj];
2572        if (iRow==otherRow) {
2573          up=false;
2574          break;
2575        }
2576      }
2577      if (up)
2578        upList[nUp++]=iColumn;
2579      else
2580        downList[nDown++]=iColumn;
2581    }
2582  }
2583  //printf("way %d\n",way);
2584  // create object
2585  //printf("would fix %d down and %d up\n",nDown,nUp);
2586  CbcBranchingObject * branch
2587     = new CbcFixingBranchingObject(model_,way,
2588                                         nDown,downList,nUp,upList);
2589  delete [] upList;
2590  delete [] downList;
2591  return branch;
2592}
2593// Default Constructor
2594CbcFixingBranchingObject::CbcFixingBranchingObject()
2595  :CbcBranchingObject()
2596{
2597  numberDown_=0;
2598  numberUp_=0;
2599  downList_=NULL;
2600  upList_=NULL;
2601}
2602
2603// Useful constructor
2604CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
2605                                                    int way ,
2606                                                    int numberOnDownSide, const int * down,
2607                                                    int numberOnUpSide, const int * up)
2608  :CbcBranchingObject(model,0,way,0.5)
2609{
2610  numberDown_=numberOnDownSide;
2611  numberUp_=numberOnUpSide;
2612  downList_ = CoinCopyOfArray(down,numberDown_);
2613  upList_ = CoinCopyOfArray(up,numberUp_);
2614}
2615
2616// Copy constructor
2617CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) :CbcBranchingObject(rhs)
2618{
2619  numberDown_=rhs.numberDown_;
2620  numberUp_=rhs.numberUp_;
2621  downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2622  upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2623}
2624
2625// Assignment operator
2626CbcFixingBranchingObject & 
2627CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject& rhs)
2628{
2629  if (this != &rhs) {
2630    CbcBranchingObject::operator=(rhs);
2631    delete [] downList_;
2632    delete [] upList_;
2633    numberDown_=rhs.numberDown_;
2634    numberUp_=rhs.numberUp_;
2635    downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2636    upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2637  }
2638  return *this;
2639}
2640CbcBranchingObject * 
2641CbcFixingBranchingObject::clone() const
2642{ 
2643  return (new CbcFixingBranchingObject(*this));
2644}
2645
2646
2647// Destructor
2648CbcFixingBranchingObject::~CbcFixingBranchingObject ()
2649{
2650  delete [] downList_;
2651  delete [] upList_;
2652}
2653double
2654CbcFixingBranchingObject::branch()
2655{
2656  decrementNumberBranchesLeft();
2657  OsiSolverInterface * solver = model_->solver();
2658  const double * columnLower = solver->getColLower();
2659  int i;
2660  // *** for way - up means fix all those in up section
2661  if (way_<0) {
2662#ifdef FULL_PRINT
2663    printf("Down Fix ");
2664#endif
2665    //printf("Down Fix %d\n",numberDown_);
2666    for (i=0;i<numberDown_;i++) {
2667      int iColumn = downList_[i];
2668      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2669#ifdef FULL_PRINT
2670      printf("Setting bound on %d to lower bound\n",iColumn);
2671#endif
2672    }
2673    way_=1;       // Swap direction
2674  } else {
2675#ifdef FULL_PRINT
2676    printf("Up Fix ");
2677#endif
2678    //printf("Up Fix %d\n",numberUp_);
2679    for (i=0;i<numberUp_;i++) {
2680      int iColumn = upList_[i];
2681      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2682#ifdef FULL_PRINT
2683      printf("Setting bound on %d to lower bound\n",iColumn);
2684#endif
2685    }
2686    way_=-1;      // Swap direction
2687  }
2688#ifdef FULL_PRINT
2689  printf("\n");
2690#endif
2691  return 0.0;
2692}
2693void
2694CbcFixingBranchingObject::print()
2695{
2696  int i;
2697  // *** for way - up means fix all those in up section
2698  if (way_<0) {
2699    printf("Down Fix ");
2700    for (i=0;i<numberDown_;i++) {
2701      int iColumn = downList_[i];
2702      printf("%d ",iColumn);
2703    }
2704  } else {
2705    printf("Up Fix ");
2706    for (i=0;i<numberUp_;i++) {
2707      int iColumn = upList_[i];
2708      printf("%d ",iColumn);
2709    }
2710  }
2711  printf("\n");
2712}
2713// Default Constructor
2714CbcNWay::CbcNWay ()
2715  : CbcObject(),
2716    numberMembers_(0),
2717    members_(NULL),
2718    consequence_(NULL)
2719{
2720}
2721
2722// Useful constructor (which are integer indices)
2723CbcNWay::CbcNWay (CbcModel * model, int numberMembers,
2724                  const int * which, int identifier)
2725  : CbcObject(model)
2726{
2727  id_=identifier;
2728  numberMembers_=numberMembers;
2729  if (numberMembers_) {
2730    members_ = new int[numberMembers_];
2731    memcpy(members_,which,numberMembers_*sizeof(int));
2732  } else {
2733    members_ = NULL;
2734  }
2735  consequence_ = NULL;
2736}
2737
2738// Copy constructor
2739CbcNWay::CbcNWay ( const CbcNWay & rhs)
2740  :CbcObject(rhs)
2741{
2742  numberMembers_ = rhs.numberMembers_;
2743  consequence_ = NULL;
2744  if (numberMembers_) {
2745    members_ = new int[numberMembers_];
2746    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
2747    if (rhs.consequence_) {
2748      consequence_ = new CbcConsequence * [numberMembers_];
2749      for (int i=0;i<numberMembers_;i++) {
2750        if (rhs.consequence_[i])
2751          consequence_[i]= rhs.consequence_[i]->clone();
2752        else
2753          consequence_[i]=NULL;
2754      }
2755    }
2756  } else {
2757    members_ = NULL;
2758  }
2759}
2760
2761// Clone
2762CbcObject *
2763CbcNWay::clone() const
2764{
2765  return new CbcNWay(*this);
2766}
2767
2768// Assignment operator
2769CbcNWay & 
2770CbcNWay::operator=( const CbcNWay& rhs)
2771{
2772  if (this!=&rhs) {
2773    CbcObject::operator=(rhs);
2774    delete [] members_;
2775    numberMembers_ = rhs.numberMembers_;
2776    if (consequence_) {
2777      for (int i=0;i<numberMembers_;i++) 
2778        delete consequence_[i];
2779      delete [] consequence_;
2780      consequence_=NULL;
2781    }
2782    if (numberMembers_) {
2783      members_ = new int[numberMembers_];
2784      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
2785    } else {
2786      members_ = NULL;
2787    }
2788    if (rhs.consequence_) {
2789      consequence_ = new CbcConsequence * [numberMembers_];
2790      for (int i=0;i<numberMembers_;i++) {
2791        if (rhs.consequence_[i])
2792          consequence_[i]= rhs.consequence_[i]->clone();
2793        else
2794          consequence_[i]=NULL;
2795      }
2796    }
2797  }
2798  return *this;
2799}
2800
2801// Destructor
2802CbcNWay::~CbcNWay ()
2803{
2804  delete [] members_;
2805  if (consequence_) {
2806    for (int i=0;i<numberMembers_;i++) 
2807      delete consequence_[i];
2808    delete [] consequence_;
2809  }
2810}
2811// Set up a consequence for a single member
2812void 
2813CbcNWay::setConsequence(int iColumn, const CbcConsequence & consequence)
2814{
2815  if (!consequence_) {
2816    consequence_ = new CbcConsequence * [numberMembers_];
2817    for (int i=0;i<numberMembers_;i++) 
2818      consequence_[i]=NULL;
2819  }
2820  for (int i=0;i<numberMembers_;i++) {
2821    if (members_[i]==iColumn) {
2822      consequence_[i]=consequence.clone();
2823      break;
2824    }
2825  }
2826}
2827
2828// Applies a consequence for a single member
2829void 
2830CbcNWay::applyConsequence(int iSequence, int state) const
2831{
2832  assert (state==-9999||state==9999);
2833  if (consequence_) {
2834    CbcConsequence * consequence = consequence_[iSequence];
2835    if (consequence) 
2836      consequence->applyToSolver(model_->solver(),state);
2837  }
2838}
2839 
2840// Infeasibility - large is 0.5
2841double 
2842CbcNWay::infeasibility(int & preferredWay) const
2843{
2844  int numberUnsatis=0;
2845  int j;
2846  OsiSolverInterface * solver = model_->solver();
2847  const double * solution = model_->testSolution();
2848  const double * lower = solver->getColLower();
2849  const double * upper = solver->getColUpper();
2850  double largestValue=0.0;
2851 
2852  double integerTolerance = 
2853    model_->getDblParam(CbcModel::CbcIntegerTolerance);
2854
2855  for (j=0;j<numberMembers_;j++) {
2856    int iColumn = members_[j];
2857    double value = solution[iColumn];
2858    value = CoinMax(value, lower[iColumn]);
2859    value = CoinMin(value, upper[iColumn]);
2860    double distance = CoinMin(value-lower[iColumn],upper[iColumn]-value);
2861    if (distance>integerTolerance) {
2862      numberUnsatis++;
2863      largestValue = CoinMax(distance,largestValue);
2864    }
2865  }
2866  preferredWay=1;
2867  if (numberUnsatis) {
2868    return largestValue;
2869  } else {
2870    return 0.0; // satisfied
2871  }
2872}
2873
2874// This looks at solution and sets bounds to contain solution
2875void 
2876CbcNWay::feasibleRegion()
2877{
2878  int j;
2879  OsiSolverInterface * solver = model_->solver();
2880  const double * solution = model_->testSolution();
2881  const double * lower = solver->getColLower();
2882  const double * upper = solver->getColUpper();
2883  double integerTolerance = 
2884    model_->getDblParam(CbcModel::CbcIntegerTolerance);
2885  for (j=0;j<numberMembers_;j++) {
2886    int iColumn = members_[j];
2887    double value = solution[iColumn];
2888    value = CoinMax(value, lower[iColumn]);
2889    value = CoinMin(value, upper[iColumn]);
2890    if (value>=upper[iColumn]-integerTolerance) {
2891      solver->setColLower(iColumn,upper[iColumn]);
2892    } else {
2893      assert (value<=lower[iColumn]+integerTolerance);
2894      solver->setColUpper(iColumn,lower[iColumn]);
2895    }
2896  }
2897}
2898// Redoes data when sequence numbers change
2899void 
2900CbcNWay::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
2901{
2902  model_=model;
2903  int n2=0;
2904  for (int j=0;j<numberMembers_;j++) {
2905    int iColumn = members_[j];
2906    int i;
2907    for (i=0;i<numberColumns;i++) {
2908      if (originalColumns[i]==iColumn)
2909        break;
2910    }
2911    if (i<numberColumns) {
2912      members_[n2]=i;
2913      consequence_[n2++]=consequence_[j];
2914    } else {
2915      delete consequence_[j];
2916    }
2917  }
2918  if (n2<numberMembers_) {
2919    printf("** NWay number of members reduced from %d to %d!\n",numberMembers_,n2);
2920    numberMembers_=n2;
2921  }
2922}
2923
2924
2925// Creates a branching object
2926CbcBranchingObject * 
2927CbcNWay::createBranch(int way) 
2928{
2929  int numberFree=0;
2930  int j;
2931
2932  OsiSolverInterface * solver = model_->solver();
2933  const double * solution = model_->testSolution();
2934  const double * lower = solver->getColLower();
2935  const double * upper = solver->getColUpper();
2936  int * list = new int[numberMembers_];
2937  double * sort = new double[numberMembers_];
2938
2939  for (j=0;j<numberMembers_;j++) {
2940    int iColumn = members_[j];
2941    double value = solution[iColumn];
2942    value = CoinMax(value, lower[iColumn]);
2943    value = CoinMin(value, upper[iColumn]);
2944    if (upper[iColumn]>lower[iColumn]) {
2945      double distance = upper[iColumn]-value;
2946      list[numberFree]=j;
2947      sort[numberFree++]=distance;
2948    }
2949  }
2950  assert (numberFree);
2951  // sort
2952  CoinSort_2(sort,sort+numberFree,list);
2953  // create object
2954  CbcBranchingObject * branch;
2955  branch = new CbcNWayBranchingObject(model_,this,numberFree,list);
2956  branch->setOriginalObject(this);
2957  delete [] list;
2958  delete [] sort;
2959  return branch;
2960}
2961 
2962// Default Constructor
2963CbcNWayBranchingObject::CbcNWayBranchingObject()
2964  :CbcBranchingObject()
2965{
2966  order_=NULL;
2967  object_=NULL;
2968  numberInSet_=0;
2969  way_=0;
2970}
2971
2972// Useful constructor
2973CbcNWayBranchingObject::CbcNWayBranchingObject (CbcModel * model,
2974                                                const CbcNWay * nway, 
2975                                                int number, const int * order)
2976  :CbcBranchingObject(model,nway->id(),-1,0.5)
2977{
2978  numberBranches_ = number;
2979  order_ = new int [number];
2980  object_=nway;
2981  numberInSet_=number;
2982  memcpy(order_,order,number*sizeof(int));
2983}
2984
2985// Copy constructor
2986CbcNWayBranchingObject::CbcNWayBranchingObject ( const CbcNWayBranchingObject & rhs) :CbcBranchingObject(rhs)
2987{
2988  numberInSet_=rhs.numberInSet_;
2989  object_=rhs.object_;
2990  if (numberInSet_) {
2991    order_ = new int [numberInSet_];
2992    memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
2993  } else {
2994    order_=NULL;
2995  }   
2996}
2997
2998// Assignment operator
2999CbcNWayBranchingObject & 
3000CbcNWayBranchingObject::operator=( const CbcNWayBranchingObject& rhs)
3001{
3002  if (this != &rhs) {
3003    CbcBranchingObject::operator=(rhs);
3004    object_=rhs.object_;
3005    delete [] order_;
3006    numberInSet_=rhs.numberInSet_;
3007    if (numberInSet_) {
3008      order_ = new int [numberInSet_];
3009      memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
3010    } else {
3011      order_=NULL;
3012    }   
3013  }
3014  return *this;
3015}
3016CbcBranchingObject * 
3017CbcNWayBranchingObject::clone() const
3018{ 
3019  return (new CbcNWayBranchingObject(*this));
3020}
3021
3022
3023// Destructor
3024CbcNWayBranchingObject::~CbcNWayBranchingObject ()
3025{
3026  delete [] order_;
3027}
3028double
3029CbcNWayBranchingObject::branch()
3030{
3031  int which = branchIndex_;
3032  branchIndex_++;
3033  assert (numberBranchesLeft()>=0);
3034  if (which==0) {
3035    // first branch so way_ may mean something
3036    assert (way_==-1||way_==1);
3037    if (way_==-1)
3038      which++;
3039  } else if (which==1) {
3040    // second branch so way_ may mean something
3041    assert (way_==-1||way_==1);
3042    if (way_==-1)
3043      which--;
3044    // switch way off
3045    way_=0;
3046  }
3047  const double * lower = model_->solver()->getColLower();
3048  const double * upper = model_->solver()->getColUpper();
3049  const int * members = object_->members();
3050  for (int j=0;j<numberInSet_;j++) {
3051    int iSequence = order_[j];
3052    int iColumn = members[iSequence];
3053    if (j!=which) {
3054      model_->solver()->setColUpper(iColumn,lower[iColumn]);
3055      //model_->solver()->setColLower(iColumn,lower[iColumn]);
3056      assert (lower[iColumn]>-1.0e20);
3057      // apply any consequences
3058      object_->applyConsequence(iSequence,-9999);
3059    } else {
3060      model_->solver()->setColLower(iColumn,upper[iColumn]);
3061      //model_->solver()->setColUpper(iColumn,upper[iColumn]);
3062#ifdef FULL_PRINT
3063      printf("Up Fix %d to %g\n",iColumn,upper[iColumn]);
3064#endif
3065      assert (upper[iColumn]<1.0e20);
3066      // apply any consequences
3067      object_->applyConsequence(iSequence,9999);
3068    }
3069  }
3070  return 0.0;
3071}
3072void
3073CbcNWayBranchingObject::print()
3074{
3075  printf("NWay - Up Fix ");
3076  const int * members = object_->members();
3077  for (int j=0;j<way_;j++) {
3078    int iColumn = members[order_[j]];
3079    printf("%d ",iColumn);
3080  }
3081  printf("\n");
3082}
3083
3084// Default Constructor
3085CbcFixVariable::CbcFixVariable ()
3086  : CbcConsequence(),
3087    numberStates_(0),
3088    states_(NULL),
3089    startLower_(NULL),
3090    startUpper_(NULL),
3091    newBound_(NULL),
3092    variable_(NULL)
3093{
3094}
3095
3096// One useful Constructor
3097CbcFixVariable::CbcFixVariable (int numberStates,const int * states, const int * numberNewLower, 
3098                                const int ** newLowerValue,
3099                                const int ** lowerColumn,
3100                                const int * numberNewUpper, const int ** newUpperValue,
3101                                const int ** upperColumn)
3102  : CbcConsequence(),
3103    states_(NULL),
3104    startLower_(NULL),
3105    startUpper_(NULL),
3106    newBound_(NULL),
3107    variable_(NULL)
3108{
3109  // How much space
3110  numberStates_ = numberStates;
3111  if (numberStates_) {
3112    states_ = new int[numberStates_];
3113    memcpy(states_,states,numberStates_*sizeof(int));
3114    int i;
3115    int n=0;
3116    startLower_ = new int[numberStates_+1];
3117    startUpper_ = new int[numberStates_+1];
3118    startLower_[0]=0;
3119    //count
3120    for (i=0;i<numberStates_;i++) {
3121      n += numberNewLower[i];
3122      startUpper_[i]=n;
3123      n += numberNewUpper[i];
3124      startLower_[i+1]=n;
3125    }
3126    newBound_ = new double [n];
3127    variable_ = new int [n];
3128    n=0;
3129    for (i=0;i<numberStates_;i++) {
3130      int j;
3131      int k;
3132      const int * bound;
3133      const int * variable;
3134      k=numberNewLower[i];
3135      bound = newLowerValue[i];
3136      variable = lowerColumn[i];
3137      for (j=0;j<k;j++) {
3138        newBound_[n]=bound[j];
3139        variable_[n++]=variable[j];
3140      }
3141      k=numberNewUpper[i];
3142      bound = newUpperValue[i];
3143      variable = upperColumn[i];
3144      for (j=0;j<k;j++) {
3145        newBound_[n]=bound[j];
3146        variable_[n++]=variable[j];
3147      }
3148    }
3149  }
3150}
3151
3152// Copy constructor
3153CbcFixVariable::CbcFixVariable ( const CbcFixVariable & rhs)
3154  :CbcConsequence(rhs)
3155{
3156  numberStates_ = rhs.numberStates_;
3157  states_ = NULL;
3158  startLower_ = NULL;
3159  startUpper_ = NULL;
3160  newBound_ = NULL;
3161  variable_ = NULL;
3162  if (numberStates_) {
3163    states_ = CoinCopyOfArray(rhs.states_,numberStates_);
3164    startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
3165    startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
3166    int n=startLower_[numberStates_];
3167    newBound_ = CoinCopyOfArray(rhs.newBound_,n);
3168    variable_ = CoinCopyOfArray(rhs.variable_,n);
3169  }
3170}
3171
3172// Clone
3173CbcConsequence *
3174CbcFixVariable::clone() const
3175{
3176  return new CbcFixVariable(*this);
3177}
3178
3179// Assignment operator
3180CbcFixVariable & 
3181CbcFixVariable::operator=( const CbcFixVariable& rhs)
3182{
3183  if (this!=&rhs) {
3184    CbcConsequence::operator=(rhs);
3185    delete [] states_;
3186    delete [] startLower_;
3187    delete [] startUpper_;
3188    delete [] newBound_;
3189    delete [] variable_;
3190    states_ = NULL;
3191    startLower_ = NULL;
3192    startUpper_ = NULL;
3193    newBound_ = NULL;
3194    variable_ = NULL;
3195    numberStates_ = rhs.numberStates_;
3196    if (numberStates_) {
3197      states_ = CoinCopyOfArray(rhs.states_,numberStates_);
3198      startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
3199      startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
3200      int n=startLower_[numberStates_];
3201      newBound_ = CoinCopyOfArray(rhs.newBound_,n);
3202      variable_ = CoinCopyOfArray(rhs.variable_,n);
3203    }
3204  }
3205  return *this;
3206}
3207
3208// Destructor
3209CbcFixVariable::~CbcFixVariable ()
3210{
3211  delete [] states_;
3212  delete [] startLower_;
3213  delete [] startUpper_;
3214  delete [] newBound_;
3215  delete [] variable_;
3216}
3217// Set up a startLower for a single member
3218void 
3219CbcFixVariable::applyToSolver(OsiSolverInterface * solver, int state) const
3220{
3221  assert (state==-9999||state==9999);
3222  // Find state
3223  int find;
3224  for (find=0;find<numberStates_;find++) 
3225    if (states_[find]==state)
3226      break;
3227  if (find==numberStates_)
3228    return;
3229  int i;
3230  // Set new lower bounds
3231  for (i=startLower_[find];i<startUpper_[find];i++) {
3232    int iColumn = variable_[i];
3233    double value = newBound_[i];
3234    double oldValue = solver->getColLower()[iColumn];
3235    //printf("for %d old lower bound %g, new %g",iColumn,oldValue,value);
3236    solver->setColLower(iColumn,CoinMax(value,oldValue));
3237    //printf(" => %g\n",solver->getColLower()[iColumn]);
3238  }
3239  // Set new upper bounds
3240  for (i=startUpper_[find];i<startLower_[find+1];i++) {
3241    int iColumn = variable_[i];
3242    double value = newBound_[i];
3243    double oldValue = solver->getColUpper()[iColumn];
3244    //printf("for %d old upper bound %g, new %g",iColumn,oldValue,value);
3245    solver->setColUpper(iColumn,CoinMin(value,oldValue));
3246    //printf(" => %g\n",solver->getColUpper()[iColumn]);
3247  }
3248}
3249
3250// Default Constructor
3251CbcDummyBranchingObject::CbcDummyBranchingObject(CbcModel * model)
3252  :CbcBranchingObject(model,0,0,0.5)
3253{
3254  setNumberBranchesLeft(1);
3255}
3256
3257
3258// Copy constructor
3259CbcDummyBranchingObject::CbcDummyBranchingObject ( const CbcDummyBranchingObject & rhs) :CbcBranchingObject(rhs)
3260{
3261}
3262
3263// Assignment operator
3264CbcDummyBranchingObject & 
3265CbcDummyBranchingObject::operator=( const CbcDummyBranchingObject& rhs)
3266{
3267  if (this != &rhs) {
3268    CbcBranchingObject::operator=(rhs);
3269  }
3270  return *this;
3271}
3272CbcBranchingObject * 
3273CbcDummyBranchingObject::clone() const
3274{ 
3275  return (new CbcDummyBranchingObject(*this));
3276}
3277
3278
3279// Destructor
3280CbcDummyBranchingObject::~CbcDummyBranchingObject ()
3281{
3282}
3283
3284/*
3285  Perform a dummy branch
3286*/
3287double
3288CbcDummyBranchingObject::branch()
3289{
3290  decrementNumberBranchesLeft();
3291  return 0.0;
3292}
3293// Print what would happen 
3294void
3295CbcDummyBranchingObject::print()
3296{
3297  printf("Dummy branch\n");
3298}
Note: See TracBrowser for help on using the repository browser.