source: trunk/Cbc/src/CbcBranchActual.cpp @ 719

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

fix longthin example

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