source: stable/2.0/Cbc/src/CbcBranchActual.cpp @ 905

Last change on this file since 905 was 905, checked in by ladanyi, 13 years ago

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

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