source: trunk/CbcBranchActual.cpp @ 295

Last change on this file since 295 was 295, checked in by forrest, 15 years ago

for ampl

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