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

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

before using osi stuff

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