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

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

many changes

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