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

Last change on this file since 523 was 523, checked in by forrest, 12 years ago

nonlinear stuff

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