source: stable/2.4/Cbc/src/CbcBranchCut.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

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