source: trunk/Cbc/src/CbcBranchCut.cpp @ 1148

Last change on this file since 1148 was 1148, checked in by forrest, 11 years ago

changes to use heuristics with SOS etc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.9 KB
Line 
1// Copyright (C) 2004, 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 "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcBranchCut.hpp"
17#include "CoinSort.hpp"
18#include "CoinError.hpp"
19
20
21/** Default Constructor
22
23*/
24CbcBranchCut::CbcBranchCut ()
25  : CbcObject()
26{
27}
28
29/* Constructor so model can be passed up
30*/ 
31CbcBranchCut::CbcBranchCut (CbcModel * model)
32  : CbcObject(model)
33{
34}
35// Copy constructor
36CbcBranchCut::CbcBranchCut ( const CbcBranchCut & rhs)
37  :CbcObject(rhs)
38
39{
40}
41
42// Clone
43CbcObject *
44CbcBranchCut::clone() const
45{
46  return new CbcBranchCut(*this);
47}
48
49// Assignment operator
50CbcBranchCut & 
51CbcBranchCut::operator=( const CbcBranchCut& rhs)
52{
53  return *this;
54}
55
56// Destructor
57CbcBranchCut::~CbcBranchCut ()
58{
59}
60
61// Infeasibility - large is 0.5
62double 
63CbcBranchCut::infeasibility(int & preferredWay) const
64{
65  throw CoinError("Use of base class","infeasibility","CbcBranchCut");
66  preferredWay=-1;
67  return 0.0;
68}
69
70// This looks at solution and sets bounds to contain solution
71/** More precisely: it first forces the variable within the existing
72    bounds, and then tightens the bounds to fix the variable at the
73    nearest integer value.
74*/
75void 
76CbcBranchCut::feasibleRegion()
77{
78}
79/* Return true if branch created by object should fix variables
80 */
81bool 
82CbcBranchCut::boundBranch() const 
83{return false;}
84
85// Creates a branching object
86CbcBranchingObject * 
87CbcBranchCut::createBranch(int way) 
88{
89  throw CoinError("Use of base class","createBranch","CbcBranchCut");
90  return new CbcCutBranchingObject();
91}
92
93
94/* Given valid solution (i.e. satisfied) and reduced costs etc
95   returns a branching object which would give a new feasible
96   point in direction reduced cost says would be cheaper.
97   If no feasible point returns null
98*/
99CbcBranchingObject * 
100CbcBranchCut::preferredNewFeasible() const
101{
102  throw CoinError("Use of base class","preferredNewFeasible","CbcBranchCut");
103  return new CbcCutBranchingObject();
104}
105 
106/* Given valid solution (i.e. satisfied) and reduced costs etc
107   returns a branching object which would give a new feasible
108   point in direction opposite to one reduced cost says would be cheaper.
109   If no feasible point returns null
110*/
111CbcBranchingObject * 
112CbcBranchCut::notPreferredNewFeasible() const 
113{
114  throw CoinError("Use of base class","notPreferredNewFeasible","CbcBranchCut");
115  return new CbcCutBranchingObject();
116}
117 
118/*
119  Bounds may be tightened, so it may be good to be able to refresh the local
120  copy of the original bounds.
121 */
122void 
123CbcBranchCut::resetBounds()
124{
125}
126
127
128// Default Constructor
129CbcCutBranchingObject::CbcCutBranchingObject()
130  :CbcBranchingObject()
131{
132  down_=OsiRowCut();
133  up_=OsiRowCut();
134  canFix_=false;
135}
136
137// Useful constructor
138CbcCutBranchingObject::CbcCutBranchingObject (CbcModel * model, 
139                                              OsiRowCut & down,
140                                              OsiRowCut &up,
141                                              bool canFix)
142  :CbcBranchingObject(model,0,-1,0.0)
143{
144  down_ = down;
145  up_ = up;
146  canFix_ = canFix;
147}
148
149// Copy constructor
150CbcCutBranchingObject::CbcCutBranchingObject ( const CbcCutBranchingObject & rhs) :CbcBranchingObject(rhs)
151{
152  down_ = rhs.down_;
153  up_ = rhs.up_;
154  canFix_ = rhs.canFix_;
155}
156
157// Assignment operator
158CbcCutBranchingObject & 
159CbcCutBranchingObject::operator=( const CbcCutBranchingObject& rhs)
160{
161  if (this != &rhs) {
162    CbcBranchingObject::operator=(rhs);
163    down_ = rhs.down_;
164    up_ = rhs.up_;
165    canFix_ = rhs.canFix_;
166  }
167  return *this;
168}
169CbcBranchingObject * 
170CbcCutBranchingObject::clone() const
171{ 
172  return (new CbcCutBranchingObject(*this));
173}
174
175
176// Destructor
177CbcCutBranchingObject::~CbcCutBranchingObject ()
178{
179}
180
181/*
182  Perform a branch by adjusting bounds and/or adding a cut. Note
183  that each arm of the branch advances the object to the next arm by
184  advancing the value of way_.
185
186  Returns change in guessed objective on next branch
187*/
188double
189CbcCutBranchingObject::branch()
190{
191  decrementNumberBranchesLeft();
192  OsiRowCut * cut;
193  if (way_<0) {
194    cut = &down_;
195    way_=1;
196  } else {
197    cut = &up_;
198    way_=-1;      // Swap direction
199  }
200  printf("CUT %s ",(way_==-1) ? "up" : "down");
201  cut->print();
202  // See if cut just fixes variables
203  double lb = cut->lb();
204  double ub = cut->ub();
205  int n=cut->row().getNumElements();
206  const int * column = cut->row().getIndices();
207  const double * element = cut->row().getElements();
208  OsiSolverInterface * solver = model_->solver();
209  const double * upper = solver->getColUpper();
210  const double * lower = solver->getColLower();
211  double low = 0.0;
212  double high=0.0;
213  for (int i=0;i<n;i++) {
214    int iColumn = column[i];
215    double value = element[i];
216    if (value>0.0) {
217      high += upper[iColumn]*value;
218      low += lower[iColumn]*value;
219    } else {
220      high += lower[iColumn]*value;
221      low += upper[iColumn]*value;
222    }
223  }
224  // leave as cut
225  //model_->setNextRowCut(*cut);
226  //return 0.0;
227  // assume cut was cunningly constructed so we need not worry too much about tolerances
228  if (low+1.0e-8>=ub&&canFix_) {
229    // fix
230    for (int i=0;i<n;i++) {
231      int iColumn = column[i];
232      double value = element[i];
233      if (value>0.0) {
234        solver->setColUpper(iColumn,lower[iColumn]);
235      } else {
236        solver->setColLower(iColumn,upper[iColumn]);
237      }
238    }
239  } else if (high-1.0e-8<=lb&&canFix_) {
240    // fix
241    for (int i=0;i<n;i++) {
242      int iColumn = column[i];
243      double value = element[i];
244      if (value>0.0) {
245        solver->setColLower(iColumn,upper[iColumn]);
246      } else {
247        solver->setColUpper(iColumn,lower[iColumn]);
248      }
249    }
250  } else {
251    // leave as cut
252    model_->setNextRowCut(*cut);
253  }
254  return 0.0;
255}
256// Print what would happen 
257void
258CbcCutBranchingObject::print()
259{
260  OsiRowCut * cut;
261  if (way_<0) {
262    cut = &down_;
263    printf("CbcCut would branch down");
264  } else {
265    cut = &up_;
266    printf("CbcCut would branch up");
267  }
268  double lb = cut->lb();
269  double ub = cut->ub();
270  int n=cut->row().getNumElements();
271  const int * column = cut->row().getIndices();
272  const double * element = cut->row().getElements();
273  if (n>5) {
274    printf(" - %d elements, lo=%g, up=%g\n",n,lb,ub);
275  } else {
276    printf(" - %g <=",lb);
277    for (int i=0;i<n;i++) {
278      int iColumn = column[i];
279      double value = element[i];
280      printf(" (%d,%g)",iColumn,value);
281    }
282    printf(" <= %g\n",ub);
283  }
284}
285
286// Return true if branch should fix variables
287bool 
288CbcCutBranchingObject::boundBranch() const
289{
290  return false;
291}
292
293/** Compare the original object of \c this with the original object of \c
294    brObj. Assumes that there is an ordering of the original objects.
295    This method should be invoked only if \c this and brObj are of the same
296    type.
297    Return negative/0/positive depending on whether \c this is
298    smaller/same/larger than the argument.
299*/
300int
301CbcCutBranchingObject::compareOriginalObject
302(const CbcBranchingObject* brObj) const
303{
304  const CbcCutBranchingObject* br =
305    dynamic_cast<const CbcCutBranchingObject*>(brObj);
306  assert(br);
307  const OsiRowCut& r0 = way_ == -1 ? down_ : up_;
308  const OsiRowCut& r1 = br->way_ == -1 ? br->down_ : br->up_;
309  return r0.row().compare(r1.row());
310}
311
312/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
313    same type and must have the same original object, but they may have
314    different feasible regions.
315    Return the appropriate CbcRangeCompare value (first argument being the
316    sub/superset if that's the case). In case of overlap (and if \c
317    replaceIfOverlap is true) replace the current branching object with one
318    whose feasible region is the overlap.
319*/
320
321CbcRangeCompare
322CbcCutBranchingObject::compareBranchingObject
323(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
324{
325  const CbcCutBranchingObject* br =
326    dynamic_cast<const CbcCutBranchingObject*>(brObj);
327  assert(br);
328  OsiRowCut& r0 = way_ == -1 ? down_ : up_;
329  const OsiRowCut& r1 = br->way_ == -1 ? br->down_ : br->up_;
330  double thisBd[2];
331  thisBd[0] = r0.lb();
332  thisBd[1] = r0.ub();
333  double otherBd[2];
334  otherBd[0] = r1.lb();
335  otherBd[1] = r1.ub();
336  CbcRangeCompare comp = CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
337  if (comp != CbcRangeOverlap || (comp==CbcRangeOverlap && !replaceIfOverlap)) {
338    return comp;
339  }
340  r0.setLb(thisBd[0]);
341  r0.setUb(thisBd[1]);
342  return comp;
343}
344
345//##############################################################################
346
347/** Default Constructor
348
349  Equivalent to an unspecified binary variable.
350*/
351CbcBranchToFixLots::CbcBranchToFixLots ()
352  : CbcBranchCut(),
353    djTolerance_(COIN_DBL_MAX),
354    fractionFixed_(1.0),
355    mark_(NULL),
356    depth_(-1),
357    numberClean_(0),
358    alwaysCreate_(false)
359{
360}
361
362/* Useful constructor - passed reduced cost tolerance and fraction we would like fixed.
363   Also depth level to do at.
364   Also passed number of 1 rows which when clean triggers fix
365   Always does if all 1 rows cleaned up and number>0 or if fraction columns reached
366   Also whether to create branch if can't reach fraction.
367*/ 
368CbcBranchToFixLots::CbcBranchToFixLots (CbcModel * model, double djTolerance,
369                                        double fractionFixed, int depth,
370                                        int numberClean,
371                                        const char * mark, bool alwaysCreate)
372  : CbcBranchCut(model)
373{
374  djTolerance_ = djTolerance;
375  fractionFixed_ = fractionFixed;
376  if (mark) {
377    int numberColumns = model->getNumCols();
378    mark_ = new char[numberColumns];
379    memcpy(mark_,mark,numberColumns);
380  } else {
381    mark_ = NULL;
382  }
383  depth_ = depth;
384  assert (model);
385  OsiSolverInterface * solver = model_->solver();
386  matrixByRow_ = *solver->getMatrixByRow();
387  numberClean_ = numberClean;
388  alwaysCreate_ = alwaysCreate;
389}
390// Copy constructor
391CbcBranchToFixLots::CbcBranchToFixLots ( const CbcBranchToFixLots & rhs)
392  :CbcBranchCut(rhs)
393{
394  djTolerance_ = rhs.djTolerance_;
395  fractionFixed_ = rhs.fractionFixed_;
396  int numberColumns = model_->getNumCols();
397  mark_ = CoinCopyOfArray(rhs.mark_,numberColumns);
398  matrixByRow_=rhs.matrixByRow_;
399  depth_ = rhs.depth_;
400  numberClean_ = rhs.numberClean_;
401  alwaysCreate_ = rhs.alwaysCreate_;
402}
403
404// Clone
405CbcObject *
406CbcBranchToFixLots::clone() const
407{
408  return new CbcBranchToFixLots(*this);
409}
410
411// Assignment operator
412CbcBranchToFixLots & 
413CbcBranchToFixLots::operator=( const CbcBranchToFixLots& rhs)
414{
415  if (this!=&rhs) {
416    CbcBranchCut::operator=(rhs);
417    djTolerance_ = rhs.djTolerance_;
418    fractionFixed_ = rhs.fractionFixed_;
419    int numberColumns = model_->getNumCols();
420    delete [] mark_;
421    mark_ = CoinCopyOfArray(rhs.mark_,numberColumns);
422    matrixByRow_=rhs.matrixByRow_;
423    depth_ = rhs.depth_;
424    numberClean_ = rhs.numberClean_;
425    alwaysCreate_ = rhs.alwaysCreate_;
426  }
427  return *this;
428}
429
430// Destructor
431CbcBranchToFixLots::~CbcBranchToFixLots ()
432{
433  delete [] mark_;
434}
435// Creates a branching object
436CbcBranchingObject * 
437CbcBranchToFixLots::createBranch(int way) 
438{
439  // by default way must be -1
440  assert (way==-1);
441  OsiSolverInterface * solver = model_->solver();
442  const double * solution = model_->testSolution();
443  const double * lower = solver->getColLower();
444  const double * upper = solver->getColUpper();
445  const double * dj = solver->getReducedCost();
446  int i;
447  int numberIntegers = model_->numberIntegers();
448  const int * integerVariable = model_->integerVariable();
449  double integerTolerance = 
450    model_->getDblParam(CbcModel::CbcIntegerTolerance);
451  // make smaller ?
452  double tolerance = CoinMin(1.0e-8,integerTolerance);
453  // How many fixed are we aiming at
454  int wantedFixed = static_cast<int> (static_cast<double>(numberIntegers)*fractionFixed_);
455  int nSort=0;
456  int numberFixed=0;
457  int numberColumns = solver->getNumCols();
458  int * sort = new int[numberColumns];
459  double * dsort = new double[numberColumns];
460  if (djTolerance_!=-1.234567) {
461    int type = shallWe();
462    assert (type);
463    // Take clean first
464    if (type==1) {
465      for (i=0;i<numberIntegers;i++) {
466        int iColumn = integerVariable[i];
467        if (upper[iColumn]>lower[iColumn]) {
468          if (!mark_||!mark_[iColumn]) {
469            if(solution[iColumn]<lower[iColumn]+tolerance) {
470              if (dj[iColumn]>djTolerance_) {
471                dsort[nSort]=-dj[iColumn];
472                sort[nSort++]=iColumn;
473              }
474            } else if (solution[iColumn]>upper[iColumn]-tolerance) {
475              if (dj[iColumn]<-djTolerance_) {
476                dsort[nSort]=dj[iColumn];
477                sort[nSort++]=iColumn;
478              }
479            }
480          }
481        } else {
482          numberFixed++;
483        }
484      }
485      // sort
486      CoinSort_2(dsort,dsort+nSort,sort);
487      nSort= CoinMin(nSort,wantedFixed-numberFixed);
488    } else if (type<10) {
489      int i;
490      //const double * rowLower = solver->getRowLower();
491      const double * rowUpper = solver->getRowUpper();
492      // Row copy
493      const double * elementByRow = matrixByRow_.getElements();
494      const int * column = matrixByRow_.getIndices();
495      const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
496      const int * rowLength = matrixByRow_.getVectorLengths();
497      const double * columnLower = solver->getColLower();
498      const double * columnUpper = solver->getColUpper();
499      const double * solution = solver->getColSolution();
500      int numberColumns = solver->getNumCols();
501      int numberRows = solver->getNumRows();
502      for (i=0;i<numberColumns;i++) {
503        sort[i]=i;
504        if (columnLower[i]!=columnUpper[i]){
505          dsort[i]=1.0e100;
506        } else {
507          dsort[i]=1.0e50;
508          numberFixed++;
509        }
510      }
511      for (i=0;i<numberRows;i++) {
512        double rhsValue = rowUpper[i];
513        bool oneRow=true;
514        // check elements
515        int numberUnsatisfied=0;
516        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
517          int iColumn = column[j];
518          double value = elementByRow[j];
519          double solValue = solution[iColumn];
520          if (columnLower[iColumn]!=columnUpper[iColumn]) {
521            if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
522              numberUnsatisfied++;
523            if (value!=1.0) {
524              oneRow=false;
525              break;
526            }
527          } else {
528            rhsValue -= value*floor(solValue+0.5);
529          }
530        }
531        if (oneRow&&rhsValue<=1.0+tolerance) {
532          if (!numberUnsatisfied) {
533            for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
534              int iColumn = column[j];
535              if (dsort[iColumn]>1.0e50){
536                dsort[iColumn]=0;
537                nSort++;
538              }
539            }
540          }
541        }
542      }
543      // sort
544      CoinSort_2(dsort,dsort+numberColumns,sort);
545    } else {
546      // new way
547      for (i=0;i<numberIntegers;i++) {
548        int iColumn = integerVariable[i];
549        if (upper[iColumn]>lower[iColumn]) {
550          if (!mark_||!mark_[iColumn]) {
551            double distanceDown=solution[iColumn]-lower[iColumn];
552            double distanceUp=upper[iColumn]-solution[iColumn];
553            double distance = CoinMin(distanceDown,distanceUp);
554            if (distance>0.001&&distance<0.5) {
555              dsort[nSort]=distance;
556              sort[nSort++]=iColumn;
557            }
558          }
559        }
560      }
561      // sort
562      CoinSort_2(dsort,dsort+nSort,sort);
563      int n=0;
564      double sum=0.0;
565      for (int k=0;k<nSort;k++) {
566        sum += dsort[k];
567        if (sum<=djTolerance_)
568          n=k;
569        else
570          break;
571      }
572      nSort = CoinMin(n,numberClean_/1000000);
573    }
574  } else {
575#define FIX_IF_LESS -0.1
576    // 3 in same row and sum <FIX_IF_LESS?
577    int numberRows = matrixByRow_.getNumRows();
578    const double * solution = model_->testSolution();
579    const int * column = matrixByRow_.getIndices();
580    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
581    const int * rowLength = matrixByRow_.getVectorLengths();
582    double bestSum=1.0;
583    int nBest=-1;
584    int kRow=-1;
585    OsiSolverInterface * solver = model_->solver();
586    for (int i=0;i<numberRows;i++) {
587      int numberUnsatisfied=0;
588      double sum=0.0;
589      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
590        int iColumn = column[j];
591        if (solver->isInteger(iColumn)) {
592          double solValue = solution[iColumn];
593          if (solValue>1.0e-5&&solValue<FIX_IF_LESS) {
594            numberUnsatisfied++;
595            sum += solValue;
596          }
597        }
598      }
599      if (numberUnsatisfied>=3&&sum<FIX_IF_LESS) {
600        // possible
601        if (numberUnsatisfied>nBest||
602            numberUnsatisfied==nBest&&sum<bestSum) {
603          nBest=numberUnsatisfied;
604          bestSum=sum;
605          kRow=i;
606        }
607      }
608    }
609    assert (nBest>0);
610    for (int j=rowStart[kRow];j<rowStart[kRow]+rowLength[kRow];j++) {
611      int iColumn = column[j];
612      if (solver->isInteger(iColumn)) {
613        double solValue = solution[iColumn];
614        if (solValue>1.0e-5&&solValue<FIX_IF_LESS) {
615          sort[nSort++]=iColumn;
616        }
617      }
618    }
619  }
620  OsiRowCut down;
621  down.setLb(-COIN_DBL_MAX);
622  double rhs=0.0;
623  for (i=0;i<nSort;i++) {
624    int iColumn = sort[i];
625    double distanceDown=solution[iColumn]-lower[iColumn];
626    double distanceUp=upper[iColumn]-solution[iColumn];
627    if(distanceDown<distanceUp) {
628      rhs += lower[iColumn];
629      dsort[i]=1.0;
630    } else {
631      rhs -= upper[iColumn];
632      dsort[i]=-1.0;
633    }
634  }
635  down.setUb(rhs);
636  down.setRow(nSort,sort,dsort);
637  down.setEffectiveness(COIN_DBL_MAX); // so will persist
638  delete [] sort;
639  delete [] dsort;
640  // up is same - just with rhs changed
641  OsiRowCut up = down;
642  up.setLb(rhs +1.0);
643  up.setUb(COIN_DBL_MAX);
644  // Say can fix one way
645  CbcCutBranchingObject * newObject = 
646    new CbcCutBranchingObject(model_,down,up,true);
647  if (model_->messageHandler()->logLevel()>1)
648    printf("creating cut in CbcBranchCut\n");
649  return newObject;
650}
651/* Does a lot of the work,
652   Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
653   10 if branching on ones away from bound
654*/
655int 
656CbcBranchToFixLots::shallWe() const
657{
658  int returnCode=0;
659  OsiSolverInterface * solver = model_->solver();
660  int numberRows = matrixByRow_.getNumRows();
661  //if (numberRows!=solver->getNumRows())
662  //return 0;
663  const double * solution = model_->testSolution();
664  const double * lower = solver->getColLower();
665  const double * upper = solver->getColUpper();
666  const double * dj = solver->getReducedCost();
667  int i;
668  int numberIntegers = model_->numberIntegers();
669  const int * integerVariable = model_->integerVariable();
670  if (numberClean_>1000000) {
671    int wanted = numberClean_%1000000;
672    int * sort = new int[numberIntegers];
673    double * dsort = new double[numberIntegers];
674    int nSort=0;
675    for (i=0;i<numberIntegers;i++) {
676      int iColumn = integerVariable[i];
677      if (upper[iColumn]>lower[iColumn]) {
678        if (!mark_||!mark_[iColumn]) {
679          double distanceDown=solution[iColumn]-lower[iColumn];
680          double distanceUp=upper[iColumn]-solution[iColumn];
681          double distance = CoinMin(distanceDown,distanceUp);
682          if (distance>0.001&&distance<0.5) {
683            dsort[nSort]=distance;
684            sort[nSort++]=iColumn;
685          }
686        }
687      }
688    }
689    // sort
690    CoinSort_2(dsort,dsort+nSort,sort);
691    int n=0;
692    double sum=0.0;
693    for (int k=0;k<nSort;k++) {
694      sum += dsort[k];
695      if (sum<=djTolerance_)
696        n=k;
697      else
698        break;
699    }
700    delete [] sort;
701    delete [] dsort;
702    return (n>=wanted) ? 10 : 0;
703  }
704  double integerTolerance = 
705    model_->getDblParam(CbcModel::CbcIntegerTolerance);
706  // make smaller ?
707  double tolerance = CoinMin(1.0e-8,integerTolerance);
708  // How many fixed are we aiming at
709  int wantedFixed = static_cast<int> (static_cast<double>(numberIntegers)*fractionFixed_);
710  if (djTolerance_<1.0e10) {
711    int nSort=0;
712    int numberFixed=0;
713    for (i=0;i<numberIntegers;i++) {
714      int iColumn = integerVariable[i];
715      if (upper[iColumn]>lower[iColumn]) {
716        if (!mark_||!mark_[iColumn]) {
717          if(solution[iColumn]<lower[iColumn]+tolerance) {
718            if (dj[iColumn]>djTolerance_) {
719              nSort++;
720            }
721          } else if (solution[iColumn]>upper[iColumn]-tolerance) {
722            if (dj[iColumn]<-djTolerance_) {
723              nSort++;
724            }
725          }
726        }
727      } else {
728        numberFixed++;
729      }
730    }
731    if (numberFixed+nSort<wantedFixed&&!alwaysCreate_) {
732      returnCode = 0;
733    } else if (numberFixed<wantedFixed) {
734      returnCode = 1;
735    } else {
736      returnCode = 0;
737    }
738  }
739  if (numberClean_) {
740    // see how many rows clean
741    int i;
742    //const double * rowLower = solver->getRowLower();
743    const double * rowUpper = solver->getRowUpper();
744    // Row copy
745    const double * elementByRow = matrixByRow_.getElements();
746    const int * column = matrixByRow_.getIndices();
747    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
748    const int * rowLength = matrixByRow_.getVectorLengths();
749    const double * columnLower = solver->getColLower();
750    const double * columnUpper = solver->getColUpper();
751    const double * solution = solver->getColSolution();
752    int numberClean=0;
753    bool someToDoYet=false;
754    int numberColumns = solver->getNumCols();
755    char * mark = new char[numberColumns];
756    int numberFixed=0;
757    for (i=0;i<numberColumns;i++) {
758      if (columnLower[i]!=columnUpper[i]){
759        mark[i]=0;
760      } else {
761        mark[i]=1;
762        numberFixed++;
763      }
764    }
765    int numberNewFixed=0;
766    for (i=0;i<numberRows;i++) {
767      double rhsValue = rowUpper[i];
768      bool oneRow=true;
769      // check elements
770      int numberUnsatisfied=0;
771      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
772        int iColumn = column[j];
773        double value = elementByRow[j];
774        double solValue = solution[iColumn];
775        if (columnLower[iColumn]!=columnUpper[iColumn]) {
776          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
777            numberUnsatisfied++;
778          if (value!=1.0) {
779            oneRow=false;
780            break;
781          }
782        } else {
783          rhsValue -= value*floor(solValue+0.5);
784        }
785      }
786      if (oneRow&&rhsValue<=1.0+tolerance) {
787        if (numberUnsatisfied) {
788          someToDoYet=true;
789        } else {
790          numberClean++;
791          for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
792            int iColumn = column[j];
793            if (columnLower[iColumn]!=columnUpper[iColumn]&&!mark[iColumn]){
794              mark[iColumn]=1;
795              numberNewFixed++;
796            }
797          }
798        }
799      }
800    }
801    delete [] mark;
802    //printf("%d clean, %d old fixed, %d new fixed\n",
803    //   numberClean,numberFixed,numberNewFixed);
804    if (someToDoYet&&numberClean<numberClean_
805        &&numberNewFixed+numberFixed<wantedFixed) {
806    } else if (numberFixed<wantedFixed) {
807      returnCode |= 2;
808    } else {
809    }
810  }
811  return returnCode;
812}
813// Infeasibility - large is 0.5
814double 
815CbcBranchToFixLots::infeasibility(int & preferredWay) const
816{
817  preferredWay=-1;
818  CbcNode * node = model_->currentNode();
819  int depth;
820  if (node)
821    depth=CoinMax(node->depth(),0);
822  else
823    return 0.0;
824  if (depth_<0) {
825    return 0.0;
826  } else if (depth_>0) {
827    if ((depth%depth_)!=0)
828      return 0.0;
829  }
830  if (djTolerance_!=-1.234567) {
831    if (!shallWe())
832      return 0.0;
833    else
834      return 1.0e20;
835  } else {
836    // See if 3 in same row and sum <FIX_IF_LESS?
837    int numberRows = matrixByRow_.getNumRows();
838    const double * solution = model_->testSolution();
839    const int * column = matrixByRow_.getIndices();
840    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
841    const int * rowLength = matrixByRow_.getVectorLengths();
842    double bestSum=1.0;
843    int nBest=-1;
844    OsiSolverInterface * solver = model_->solver();
845    for (int i=0;i<numberRows;i++) {
846      int numberUnsatisfied=0;
847      double sum=0.0;
848      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
849        int iColumn = column[j];
850        if (solver->isInteger(iColumn)) {
851          double solValue = solution[iColumn];
852          if (solValue>1.0e-5&&solValue<FIX_IF_LESS) {
853            numberUnsatisfied++;
854            sum += solValue;
855          }
856        }
857      }
858      if (numberUnsatisfied>=3&&sum<FIX_IF_LESS) {
859        // possible
860        if (numberUnsatisfied>nBest||
861            numberUnsatisfied==nBest&&sum<bestSum) {
862          nBest=numberUnsatisfied;
863          bestSum=sum;
864        }
865      }
866    }
867    if (nBest>0)
868      return 1.0e20;
869    else
870      return 0.0;
871  }
872}
873// Redoes data when sequence numbers change
874void 
875CbcBranchToFixLots::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
876{
877  model_=model;
878  if (mark_) {
879    OsiSolverInterface * solver = model_->solver();
880    int numberColumnsNow = solver->getNumCols();
881    char * temp = new char[numberColumnsNow];
882    memset(temp,0,numberColumnsNow);
883    for (int i=0;i<numberColumns;i++) {
884      int j = originalColumns[i];
885      temp[i]=mark_[j];
886    }
887    delete [] mark_;
888    mark_=temp;
889  }
890  OsiSolverInterface * solver = model_->solver();
891  matrixByRow_ = *solver->getMatrixByRow();
892}
893
894/** Default Constructor
895*/
896CbcBranchAllDifferent::CbcBranchAllDifferent ()
897  : CbcBranchCut(),
898    numberInSet_(0),
899    which_(NULL)
900{
901}
902
903/* Useful constructor - passed set of variables
904*/ 
905CbcBranchAllDifferent::CbcBranchAllDifferent (CbcModel * model, int numberInSet,
906                                              const int * members)
907  : CbcBranchCut(model)
908{
909  numberInSet_=numberInSet;
910  which_ = CoinCopyOfArray(members,numberInSet_);
911}
912// Copy constructor
913CbcBranchAllDifferent::CbcBranchAllDifferent ( const CbcBranchAllDifferent & rhs)
914  :CbcBranchCut(rhs)
915{
916  numberInSet_=rhs.numberInSet_;
917  which_ = CoinCopyOfArray(rhs.which_,numberInSet_);
918}
919
920// Clone
921CbcObject *
922CbcBranchAllDifferent::clone() const
923{
924  return new CbcBranchAllDifferent(*this);
925}
926
927// Assignment operator
928CbcBranchAllDifferent & 
929CbcBranchAllDifferent::operator=( const CbcBranchAllDifferent& rhs)
930{
931  if (this!=&rhs) {
932    CbcBranchCut::operator=(rhs);
933    delete [] which_;
934    numberInSet_=rhs.numberInSet_;
935    which_ = CoinCopyOfArray(rhs.which_,numberInSet_);
936  }
937  return *this;
938}
939
940// Destructor
941CbcBranchAllDifferent::~CbcBranchAllDifferent ()
942{
943  delete [] which_;
944}
945// Creates a branching object
946CbcBranchingObject * 
947CbcBranchAllDifferent::createBranch(int way) 
948{
949  // by default way must be -1
950  assert (way==-1);
951  const double * solution = model_->testSolution();
952  double * values = new double[numberInSet_];
953  int * which = new int[numberInSet_];
954  int i;
955  for (i=0;i<numberInSet_;i++) {
956    int iColumn = which_[i];
957    values[i]=solution[iColumn];
958    which[i]=iColumn;
959  }
960  CoinSort_2(values,values+numberInSet_,which);
961  double last = -1.0;
962  double closest=1.0;
963  int worst=-1;
964  for (i=0;i<numberInSet_;i++) {
965    if (values[i]-last<closest) {
966      closest=values[i]-last;
967      worst=i-1;
968    }
969    last=values[i];
970  }
971  assert (closest<=0.99999);
972  OsiRowCut down;
973  down.setLb(-COIN_DBL_MAX);
974  down.setUb(-1.0);
975  int pair[2];
976  double elements[]={1.0,-1.0};
977  pair[0]=which[worst];
978  pair[1]=which[worst+1];
979  delete [] values;
980  delete [] which;
981  down.setRow(2,pair,elements);
982  // up is same - just with rhs changed
983  OsiRowCut up = down;
984  up.setLb(1.0);
985  up.setUb(COIN_DBL_MAX);
986  // Say is not a fix type branch
987  CbcCutBranchingObject * newObject = 
988    new CbcCutBranchingObject(model_,down,up,false);
989  if (model_->messageHandler()->logLevel()>1)
990    printf("creating cut in CbcBranchCut\n");
991  return newObject;
992}
993// Infeasibility - large is 0.5
994double 
995CbcBranchAllDifferent::infeasibility(int & preferredWay) const
996{
997  preferredWay=-1;
998  //OsiSolverInterface * solver = model_->solver();
999  const double * solution = model_->testSolution();
1000  //const double * lower = solver->getColLower();
1001  //const double * upper = solver->getColUpper();
1002  double * values = new double[numberInSet_];
1003  int i;
1004  for (i=0;i<numberInSet_;i++) {
1005    int iColumn = which_[i];
1006    values[i]=solution[iColumn];
1007  }
1008  std::sort(values,values+numberInSet_);
1009  double last = -1.0;
1010  double closest=1.0;
1011  for (i=0;i<numberInSet_;i++) {
1012    if (values[i]-last<closest) {
1013      closest=values[i]-last;
1014    }
1015    last=values[i];
1016  }
1017  delete [] values;
1018  if (closest>0.99999)
1019    return 0.0;
1020  else
1021    return 0.5*(1.0-closest);
1022}
Note: See TracBrowser for help on using the repository browser.