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

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

towards common use with other solvers

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