source: trunk/Cbc/src/CbcHeuristicDive.cpp @ 944

Last change on this file since 944 was 944, checked in by jpgoncal, 12 years ago

Added reduced cost fixing.

File size: 24.1 KB
RevLine 
[907]1// Copyright (C) 2008, 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
8#include "CbcHeuristicDive.hpp"
9#include "CbcStrategy.hpp"
10#include  "CoinTime.hpp"
11
[944]12#ifdef COIN_HAS_CLP
13#include "OsiClpSolverInterface.hpp"
14#endif
15
[916]16//#define DIVE_FIX_BINARY_VARIABLES
[944]17//#define DIVE_DEBUG
[916]18
[907]19// Default Constructor
20CbcHeuristicDive::CbcHeuristicDive() 
21  :CbcHeuristic()
22{
23  // matrix and row copy will automatically be empty
24  downLocks_ =NULL;
25  upLocks_ = NULL;
26  percentageToFix_ = 0.2;
27  maxIterations_ = 100;
[944]28  maxTime_ = 600;
[907]29}
30
31// Constructor from model
32CbcHeuristicDive::CbcHeuristicDive(CbcModel & model)
33  :CbcHeuristic(model)
34{
35  downLocks_ =NULL;
36  upLocks_ = NULL;
37  // Get a copy of original matrix
38  assert(model.solver());
39  // model may have empty matrix - wait until setModel
40  const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
41  if (matrix) {
42    matrix_ = *matrix;
[916]43    matrixByRow_ = *model.solver()->getMatrixByRow();
[907]44    validate();
45  }
46  percentageToFix_ = 0.2;
47  maxIterations_ = 100;
[944]48  maxTime_ = 600;
[907]49}
50
51// Destructor
52CbcHeuristicDive::~CbcHeuristicDive ()
53{
54  delete [] downLocks_;
55  delete [] upLocks_;
56}
57
58// Create C++ lines to get to current state
59void 
60CbcHeuristicDive::generateCpp( FILE * fp, const char * heuristic) 
61{
62  // hard coded as CbcHeuristic virtual
63  CbcHeuristic::generateCpp(fp,heuristic);
64  if (percentageToFix_!=0.2)
65    fprintf(fp,"3  %s.setPercentageToFix(%.f);\n",heuristic,percentageToFix_);
66  else
67    fprintf(fp,"4  %s.setPercentageToFix(%.f);\n",heuristic,percentageToFix_);
68  if (maxIterations_!=100)
69    fprintf(fp,"3  %s.setMaxIterations(%d);\n",heuristic,maxIterations_);
70  else
71    fprintf(fp,"4  %s.setMaxIterations(%d);\n",heuristic,maxIterations_);
[944]72  if (maxTime_!=600)
[907]73    fprintf(fp,"3  %s.setMaxTime(%.2f);\n",heuristic,maxTime_);
74  else
75    fprintf(fp,"4  %s.setMaxTime(%.2f);\n",heuristic,maxTime_);
76}
77
78// Copy constructor
79CbcHeuristicDive::CbcHeuristicDive(const CbcHeuristicDive & rhs)
80:
81  CbcHeuristic(rhs),
82  matrix_(rhs.matrix_),
[916]83  matrixByRow_(rhs.matrixByRow_),
[907]84  percentageToFix_(rhs.percentageToFix_),
85  maxIterations_(rhs.maxIterations_),
86  maxTime_(rhs.maxTime_)
87{
88  if (rhs.downLocks_) {
89    int numberIntegers = model_->numberIntegers();
90    downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
91    upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
92  } else {
93    downLocks_ = NULL;
94    upLocks_ = NULL;
95  }
96}
97
98// Assignment operator
99CbcHeuristicDive & 
100CbcHeuristicDive::operator=( const CbcHeuristicDive& rhs)
101{
102  if (this!=&rhs) {
103    CbcHeuristic::operator=(rhs);
104    matrix_ = rhs.matrix_;
[916]105    matrixByRow_ = rhs.matrixByRow_;
[907]106    percentageToFix_ = rhs.percentageToFix_;
107    maxIterations_ = rhs.maxIterations_;
108    maxTime_ = rhs.maxTime_;
109    delete [] downLocks_;
110    delete [] upLocks_;
111    if (rhs.downLocks_) {
112      int numberIntegers = model_->numberIntegers();
113      downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
114      upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
115    } else {
116      downLocks_ = NULL;
117      upLocks_ = NULL;
118    }
119  }
120  return *this;
121}
122
123// Resets stuff if model changes
124void 
125CbcHeuristicDive::resetModel(CbcModel * model)
126{
127  model_=model;
128  assert(model_->solver());
129  // Get a copy of original matrix
130  const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
131  // model may have empty matrix - wait until setModel
132  if (matrix) {
133    matrix_ = *matrix;
[916]134    matrixByRow_ = *model->solver()->getMatrixByRow();
[907]135    validate();
136  }
137}
138
139// update model
140void CbcHeuristicDive::setModel(CbcModel * model)
141{
142  model_ = model;
143  assert(model_->solver());
144  // Get a copy of original matrix
145  const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
146  if (matrix) {
147    matrix_ = *matrix;
[916]148    matrixByRow_ = *model->solver()->getMatrixByRow();
[907]149    // make sure model okay for heuristic
150    validate();
151  }
152}
153
154bool CbcHeuristicDive::canHeuristicRun()
155{
156  return shouldHeurRun_randomChoice();
157}
158
[916]159struct PseudoReducedCost {
160  int var;
161  double pseudoRedCost;
162};
[907]163
[916]164inline bool compareBinaryVars(const PseudoReducedCost obj1,
165                              const PseudoReducedCost obj2)
166{
167  return obj1.pseudoRedCost > obj2.pseudoRedCost;
168}
169
[907]170// See if dive fractional will give better solution
171// Sets value of solution
172// Returns 1 if solution, 0 if not
173int
174CbcHeuristicDive::solution(double & solutionValue,
175                           double * betterSolution)
176{
[944]177#ifdef DIVE_DEBUG
178  std::cout<<"solutionValue = "<<solutionValue<<std::endl;
179#endif
[907]180  ++numCouldRun_;
181
182  // test if the heuristic can run
183  if(!canHeuristicRun())
184    return 0;
185
186#if 0
187  // See if to do
188  if (!when()||(when()%10==1&&model_->phase()!=1)||
189      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
190    return 0; // switched off
191#endif
192
[944]193#ifdef DIVE_DEBUG
194  int nRoundInfeasible = 0;
195  int nRoundFeasible = 0;
196  int reasonToStop = 0;
197#endif
198
[907]199  double time1 = CoinCpuTime();
200
201  OsiSolverInterface * solver = model_->solver()->clone();
202  const double * lower = solver->getColLower();
203  const double * upper = solver->getColUpper();
204  const double * rowLower = solver->getRowLower();
205  const double * rowUpper = solver->getRowUpper();
206  const double * solution = solver->getColSolution();
207  const double * objective = solver->getObjCoefficients();
208  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
209  double primalTolerance;
210  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
211
212  int numberRows = matrix_.getNumRows();
213  assert (numberRows<=solver->getNumRows());
214  int numberIntegers = model_->numberIntegers();
215  const int * integerVariable = model_->integerVariable();
216  double direction = solver->getObjSense(); // 1 for min, -1 for max
217  double newSolutionValue = direction*solver->getObjValue();
218  int returnCode = 0;
219  // Column copy
220  const double * element = matrix_.getElements();
221  const int * row = matrix_.getIndices();
222  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
223  const int * columnLength = matrix_.getVectorLengths();
[916]224  // Row copy
225  const double * elementByRow = matrixByRow_.getElements();
226  const int * column = matrixByRow_.getIndices();
227  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
228  const int * rowLength = matrixByRow_.getVectorLengths();
[907]229
230  // Get solution array for heuristic solution
231  int numberColumns = solver->getNumCols();
232  double * newSolution = new double [numberColumns];
233  memcpy(newSolution,solution,numberColumns*sizeof(double));
234
235  // vectors to store the latest variables fixed at their bounds
236  int* columnFixed = new int [numberIntegers];
237  double* originalBound = new double [numberIntegers];
238  bool * fixedAtLowerBound = new bool [numberIntegers];
[916]239  PseudoReducedCost * candidate = NULL;
240  if(binVarIndex_.size())
241    candidate = new PseudoReducedCost [binVarIndex_.size()];
[907]242
243  const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
244
245  // count how many fractional variables
246  int numberFractionalVariables = 0;
247  for (int i=0; i<numberIntegers; i++) {
248    int iColumn = integerVariable[i];
249    double value=newSolution[iColumn];
250    if (fabs(floor(value+0.5)-value)>integerTolerance) {
251      numberFractionalVariables++;
252    }
253  }
254
[916]255  const double* reducedCost = solver->getReducedCost();
256
[907]257  int iteration = 0;
258  while(numberFractionalVariables) {
259    iteration++;
260
261    // select a fractional variable to bound
[944]262    int bestColumn = -1;
[907]263    int bestRound; // -1 rounds down, +1 rounds up
[944]264    bool canRound = selectVariableToBranch(solver, newSolution, 
265                                           bestColumn, bestRound);
[907]266
[944]267    // if the solution is not trivially roundable, we don't try to round;
268    // if the solution is trivially roundable, we try to round. However,
269    // if the rounded solution is worse than the current incumbent,
270    // then we don't round and proceed normally. In this case, the
271    // bestColumn will be a trivially roundable variable
272    if(canRound) {
273      // check if by rounding all fractional variables
274      // we get a solution with an objective value
275      // better than the current best integer solution
276      double delta = 0.0;
277      for (int i=0; i<numberIntegers; i++) {
278        int iColumn = integerVariable[i];
279        double value=newSolution[iColumn];
280        if (fabs(floor(value+0.5)-value)>integerTolerance) {
281          assert(downLocks_[i]==0 || upLocks_[i]==0);
282          double obj = objective[iColumn];
283          if(downLocks_[i]==0 && upLocks_[i]==0) {
284            if(direction * obj >= 0.0)
285              delta += (floor(value) - value) * obj;
286            else
287              delta += (ceil(value) - value) * obj;
288          }
289          else if(downLocks_[i]==0)
290            delta += (floor(value) - value) * obj;
291          else
292            delta += (ceil(value) - value) * obj;
293        }
294      }
295      if(direction*(solver->getObjValue()+delta) < solutionValue) {
296#ifdef DIVE_DEBUG
297        nRoundFeasible++;
298#endif
299        // Round all the fractional variables
300        for (int i=0; i<numberIntegers; i++) {
301          int iColumn = integerVariable[i];
302          double value=newSolution[iColumn];
303          if (fabs(floor(value+0.5)-value)>integerTolerance) {
304            assert(downLocks_[i]==0 || upLocks_[i]==0);
305            if(downLocks_[i]==0 && upLocks_[i]==0) {
306              if(direction * objective[iColumn] >= 0.0)
307                newSolution[iColumn] = floor(value);
308              else
309                newSolution[iColumn] = ceil(value);
310            }
311            else if(downLocks_[i]==0)
312              newSolution[iColumn] = floor(value);
313            else
314              newSolution[iColumn] = ceil(value);
315          }
316        }
317        break;
318      }
319#ifdef DIVE_DEBUG
320      else
321        nRoundInfeasible++;
322#endif
323    }
324
325    // do reduced cost fixing
326    int numberFixed = reducedCostFix(solver);
327#ifdef DIVE_DEBUG
328    std::cout<<"numberReducedCostFixed = "<<numberFixed<<std::endl;
329#endif
[907]330     
[944]331    int numberAtBoundFixed = 0;
332#ifdef DIVE_FIX_BINARY_VARIABLES
[916]333    // fix binary variables based on pseudo reduced cost
334    if(binVarIndex_.size()) {
335      int cnt = 0;
336      for (int j=0; j<(int)binVarIndex_.size(); j++) {
[917]337        int iColumn1 = binVarIndex_[j];
338        double value = newSolution[iColumn1];
339        double maxPseudoReducedCost = 0.0;
340        if(fabs(value)<=integerTolerance &&
341           lower[iColumn1] != upper[iColumn1]) {
[944]342#ifdef DIVE_DEBUG
343          std::cout<<"iColumn1 = "<<iColumn1<<", value = "<<value<<std::endl;
344#endif
[917]345          int iRow = vbRowIndex_[j];
346          for (int k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {
347            int iColumn2 = column[k];
[944]348#ifdef DIVE_DEBUG
349            std::cout<<"iColumn2 = "<<iColumn2<<std::endl;
350#endif
[917]351            if(iColumn1 != iColumn2) {
352              double pseudoReducedCost = fabs(reducedCost[iColumn2] *
353                                              elementByRow[iColumn2] / 
354                                              elementByRow[iColumn1]);
[944]355#ifdef DIVE_DEBUG
356              std::cout<<"reducedCost["<<iColumn2<<"] = "
357                       <<reducedCost[iColumn2]
358                       <<", elementByRow["<<iColumn2<<"] = "<<elementByRow[iColumn2]
359                       <<", elementByRow["<<iColumn1<<"] = "<<elementByRow[iColumn1]
360                       <<", pseudoRedCost = "<<pseudoReducedCost
361                       <<std::endl;
362#endif
[917]363              if(pseudoReducedCost > maxPseudoReducedCost)
364                maxPseudoReducedCost = pseudoReducedCost;
365            }
366          }
[944]367#ifdef DIVE_DEBUG
368          std::cout<<", maxPseudoRedCost = "<<maxPseudoReducedCost<<std::endl;
369#endif
[917]370          candidate[cnt].var = iColumn1;
371          candidate[cnt++].pseudoRedCost = maxPseudoReducedCost;
372        }
373      }
[944]374#ifdef DIVE_DEBUG
375      std::cout<<"candidates for rounding = "<<cnt<<std::endl;
376#endif
[917]377      std::sort(candidate, candidate+cnt, compareBinaryVars);
378      for (int i=0; i<cnt; i++) {
379        int iColumn = candidate[i].var;
380        if (numberAtBoundFixed < maxNumberAtBoundToFix) {
381          columnFixed[numberAtBoundFixed] = iColumn;
382          originalBound[numberAtBoundFixed] = upper[iColumn];
383          fixedAtLowerBound[numberAtBoundFixed] = true;
384          solver->setColUpper(iColumn, lower[iColumn]);
385          numberAtBoundFixed++;
386          if(numberAtBoundFixed == maxNumberAtBoundToFix)
387            break;
388        }
389      }
390    }
391#endif
[916]392
[944]393    // fix other integer variables that are at their bounds
[907]394    for (int i=0; i<numberIntegers; i++) {
395      int iColumn = integerVariable[i];
396      double value=newSolution[iColumn];
397      if(fabs(floor(value+0.5)-value)<=integerTolerance && 
398         numberAtBoundFixed < maxNumberAtBoundToFix) {
399        // fix the variable at one of its bounds
400        if (fabs(lower[iColumn]-value)<=integerTolerance &&
401            lower[iColumn] != upper[iColumn]) {
402          columnFixed[numberAtBoundFixed] = iColumn;
403          originalBound[numberAtBoundFixed] = upper[iColumn];
404          fixedAtLowerBound[numberAtBoundFixed] = true;
405          solver->setColUpper(iColumn, lower[iColumn]);
406          numberAtBoundFixed++;
407        }
408        else if(fabs(upper[iColumn]-value)<=integerTolerance &&
409            lower[iColumn] != upper[iColumn]) {
410          columnFixed[numberAtBoundFixed] = iColumn;
411          originalBound[numberAtBoundFixed] = lower[iColumn];
412          fixedAtLowerBound[numberAtBoundFixed] = false;
413          solver->setColLower(iColumn, upper[iColumn]);
414          numberAtBoundFixed++;
415        }
416        if(numberAtBoundFixed == maxNumberAtBoundToFix)
417          break;
418      }
419    }
[944]420#ifdef DIVE_DEBUG
421    std::cout<<"numberAtBoundFixed = "<<numberAtBoundFixed<<std::endl;
422#endif
[907]423
424    double originalBoundBestColumn;
425    if(bestColumn >= 0) {
426      if(bestRound < 0) {
427        originalBoundBestColumn = upper[bestColumn];
428        solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
429      }
430      else {
431        originalBoundBestColumn = lower[bestColumn];
432        solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
433      }
434    } else {
435      break;
436    }
437    int originalBestRound = bestRound;
438    while (1) {
439
440      solver->resolve();
441
442      if(!solver->isProvenOptimal()) {
443        if(numberAtBoundFixed > 0) {
444          // Remove the bound fix for variables that were at bounds
445          for(int i=0; i<numberAtBoundFixed; i++) {
446            int iColFixed = columnFixed[i];
447            if(fixedAtLowerBound[i])
448              solver->setColUpper(iColFixed, originalBound[i]);
449            else
450              solver->setColLower(iColFixed, originalBound[i]);
451          }
452          numberAtBoundFixed = 0;
453        }
454        else if(bestRound == originalBestRound) {
455          bestRound *= (-1);
456          if(bestRound < 0) {
457            solver->setColLower(bestColumn, originalBoundBestColumn);
458            solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
459          }
460          else {
461            solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
462            solver->setColUpper(bestColumn, originalBoundBestColumn);
463          }
464        }
465        else
466          break;
467      }
468      else
469        break;
470    }
471
[920]472    if(!solver->isProvenOptimal() || 
[944]473       direction*solver->getObjValue() >= solutionValue) {
474#ifdef DIVE_DEBUG
475      reasonToStop = 1;
476#endif
[907]477      break;
[944]478    }
[907]479
480    if(iteration > maxIterations_) {
[944]481#ifdef DIVE_DEBUG
482      reasonToStop = 2;
483#endif
[907]484      break;
485    }
486
487    if(CoinCpuTime()-time1 > maxTime_) {
[944]488#ifdef DIVE_DEBUG
489      reasonToStop = 3;
490#endif
[907]491      break;
492    }
493
494    memcpy(newSolution,solution,numberColumns*sizeof(double));
495    numberFractionalVariables = 0;
496    for (int i=0; i<numberIntegers; i++) {
497      int iColumn = integerVariable[i];
498      double value=newSolution[iColumn];
499      if (fabs(floor(value+0.5)-value)>integerTolerance) {
500        numberFractionalVariables++;
501      }
502    }
503
504  }
505
506
507  double * rowActivity = new double[numberRows];
508  memset(rowActivity,0,numberRows*sizeof(double));
509
510  // re-compute new solution value
511  double objOffset=0.0;
512  solver->getDblParam(OsiObjOffset,objOffset);
513  newSolutionValue = -objOffset;
514  for (int i=0 ; i<numberColumns ; i++ )
515    newSolutionValue += objective[i]*newSolution[i];
516  newSolutionValue *= direction;
517    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
518  if (newSolutionValue<solutionValue) {
519    // paranoid check
520    memset(rowActivity,0,numberRows*sizeof(double));
521    for (int i=0;i<numberColumns;i++) {
522      int j;
523      double value = newSolution[i];
524      if (value) {
525        for (j=columnStart[i];
526             j<columnStart[i]+columnLength[i];j++) {
527          int iRow=row[j];
528          rowActivity[iRow] += value*element[j];
529        }
530      }
531    }
532    // check was approximately feasible
533    bool feasible=true;
534    for (int i=0;i<numberRows;i++) {
535      if(rowActivity[i]<rowLower[i]) {
536        if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
537          feasible = false;
538      } else if(rowActivity[i]>rowUpper[i]) {
539        if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
540          feasible = false;
541      }
542    }
543    for (int i=0; i<numberIntegers; i++) {
544      int iColumn = integerVariable[i];
545      double value=newSolution[iColumn];
546      if (fabs(floor(value+0.5)-value)>integerTolerance) {
547        feasible = false;
548        break;
549      }
550    }
551    if (feasible) {
552      // new solution
553      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
554      solutionValue = newSolutionValue;
555      //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
556      returnCode=1;
557    } else {
558      // Can easily happen
559      //printf("Debug CbcHeuristicDive giving bad solution\n");
560    }
561  }
562
[944]563#ifdef DIVE_DEBUG
564  std::cout<<"nRoundInfeasible = "<<nRoundInfeasible
565           <<", nRoundFeasible = "<<nRoundFeasible
566           <<", returnCode = "<<returnCode
567           <<", reasonToStop = "<<reasonToStop
568           <<", iterations = "<<iteration<<std::endl;
569#endif
570
[907]571  delete [] newSolution;
572  delete [] columnFixed;
573  delete [] originalBound;
574  delete [] fixedAtLowerBound;
[916]575  delete [] candidate;
[907]576  delete [] rowActivity;
577  delete solver;
578  return returnCode;
579}
580
581// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
582void 
583CbcHeuristicDive::validate() 
584{
585  if (model_&&when()<10) {
586    if (model_->numberIntegers()!=
587        model_->numberObjects())
588      setWhen(0);
589  }
590
591  int numberIntegers = model_->numberIntegers();
592  const int * integerVariable = model_->integerVariable();
593  delete [] downLocks_;
594  delete [] upLocks_;
595  downLocks_ = new unsigned short [numberIntegers];
596  upLocks_ = new unsigned short [numberIntegers];
597  // Column copy
598  const double * element = matrix_.getElements();
599  const int * row = matrix_.getIndices();
600  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
601  const int * columnLength = matrix_.getVectorLengths();
602  const double * rowLower = model_->solver()->getRowLower();
603  const double * rowUpper = model_->solver()->getRowUpper();
604  for (int i=0;i<numberIntegers;i++) {
605    int iColumn = integerVariable[i];
606    int down=0;
607    int up=0;
608    if (columnLength[iColumn]>65535) {
609      setWhen(0);
610      break; // unlikely to work
611    }
612    for (CoinBigIndex j=columnStart[iColumn];
613         j<columnStart[iColumn]+columnLength[iColumn];j++) {
614      int iRow=row[j];
615      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
616        up++;
617        down++;
618      } else if (element[j]>0.0) {
619        if (rowUpper[iRow]<1.0e20)
620          up++;
621        else
622          down++;
623      } else {
624        if (rowLower[iRow]>-1.0e20)
625          up++;
626        else
627          down++;
628      }
629    }
630    downLocks_[i] = (unsigned short) down;
631    upLocks_[i] = (unsigned short) up;
632  }
[916]633
634#ifdef DIVE_FIX_BINARY_VARIABLES
635  selectBinaryVariables();
636#endif
[907]637}
[916]638
639// Select candidate binary variables for fixing
640void
641CbcHeuristicDive::selectBinaryVariables()
642{
643  // Row copy
644  const double * elementByRow = matrixByRow_.getElements();
645  const int * column = matrixByRow_.getIndices();
646  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
647  const int * rowLength = matrixByRow_.getVectorLengths();
648
649  const int numberRows = matrixByRow_.getNumRows();
650  const int numberCols = matrixByRow_.getNumCols();
651
652  const double * lower = model_->solver()->getColLower();
653  const double * upper = model_->solver()->getColUpper();
654  const double * rowLower = model_->solver()->getRowLower();
655  const double * rowUpper = model_->solver()->getRowUpper();
656
657  //  const char * integerType = model_->integerType();
658 
659
660  //  const int numberIntegers = model_->numberIntegers();
661  //  const int * integerVariable = model_->integerVariable();
662  const double * objective = model_->solver()->getObjCoefficients();
663
664  // vector to store the row number of variable bound rows
665  int* rowIndexes = new int [numberCols];
666  memset(rowIndexes, -1, numberCols*sizeof(int));
667
668  for(int i=0; i<numberRows; i++) {
[944]669    int positiveBinary = -1;
670    int negativeBinary = -1;
671    int nPositiveOther = 0;
672    int nNegativeOther = 0;
[917]673    for (int k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
674      int iColumn = column[k];
[944]675      if(model_->solver()->isInteger(iColumn) &&
676         lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
677         objective[iColumn] == 0.0 &&
678         elementByRow[k] > 0.0 &&
679         positiveBinary < 0)
680        positiveBinary = iColumn;
681      else if(model_->solver()->isInteger(iColumn) &&
682              lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
683              objective[iColumn] == 0.0 &&
684              elementByRow[k] < 0.0 &&
685              negativeBinary < 0)
686        negativeBinary = iColumn;
687      else if((elementByRow[k]> 0.0 &&
688               lower[iColumn] >= 0.0) ||
689              (elementByRow[k] < 0.0 &&
690               upper[iColumn] <= 0.0))
691        nPositiveOther++;
692      else if((elementByRow[k]> 0.0 &&
693               lower[iColumn] <= 0.0) ||
694              (elementByRow[k] < 0.0 &&
695               upper[iColumn] >= 0.0))
696        nNegativeOther++;
697      if(nPositiveOther > 0 && nNegativeOther > 0)
698        break;
[917]699    }
[944]700    int binVar = -1;
701    if(positiveBinary >= 0 &&
702       (negativeBinary >= 0 || nNegativeOther > 0) &&
703       nPositiveOther == 0 &&
704       rowLower[i] == 0.0 &&
705       rowUpper[i] > 0.0)
706      binVar = positiveBinary;
707    else if(negativeBinary >= 0 &&
708            (positiveBinary >= 0 || nPositiveOther > 0) &&
709            nNegativeOther == 0 &&
710            rowLower[i] < 0.0 &&
711            rowUpper[i] == 0.0)
712      binVar = negativeBinary;
713    if(binVar >= 0) {
[917]714      if(rowIndexes[binVar] == -1)
715        rowIndexes[binVar] = i;
716      else if(rowIndexes[binVar] >= 0)
717        rowIndexes[binVar] = -2;
718    }
719  }
720
721  for(int j=0; j<numberCols; j++) {
722    if(rowIndexes[j] >= 0) {
723      binVarIndex_.push_back(j);
724      vbRowIndex_.push_back(rowIndexes[j]);
725    }
726  }
727
[944]728#ifdef DIVE_DEBUG
[917]729  std::cout<<"number vub Binary = "<<binVarIndex_.size()<<std::endl;
[944]730#endif
[917]731
732  delete [] rowIndexes;
733   
734}
735
[944]736/*
737  Perform reduced cost fixing on integer variables.
[917]738
[944]739  The variables in question are already nonbasic at bound. We're just nailing
740  down the current situation.
741*/
[917]742
[944]743int CbcHeuristicDive::reducedCostFix (OsiSolverInterface* solver)
[917]744
[944]745{
746#if 0
747  if(!solverCharacteristics_->reducedCostsAccurate())
748    return 0; //NLP
749#endif
750  double cutoff = model_->getCutoff() ;
751#ifdef DIVE_DEBUG
752  std::cout<<"cutoff = "<<cutoff<<std::endl;
753#endif
754  double direction = solver->getObjSense() ;
755  double gap = cutoff - solver->getObjValue()*direction ;
756  double tolerance;
757  solver->getDblParam(OsiDualTolerance,tolerance) ;
758  if (gap<=0.0)
759    gap = tolerance; //return 0;
760  gap += 100.0*tolerance;
761  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
[917]762
[944]763  const double *lower = solver->getColLower() ;
764  const double *upper = solver->getColUpper() ;
765  const double *solution = solver->getColSolution() ;
766  const double *reducedCost = solver->getReducedCost() ;
[917]767
[944]768  int numberIntegers = model_->numberIntegers();
769  const int * integerVariable = model_->integerVariable();
[917]770
[944]771  int numberFixed = 0 ;
[917]772
[944]773# ifdef COIN_HAS_CLP
774  OsiClpSolverInterface * clpSolver
775    = dynamic_cast<OsiClpSolverInterface *> (solver);
776  ClpSimplex * clpSimplex=NULL;
777  if (clpSolver) 
778    clpSimplex = clpSolver->getModelPtr();
779# endif
780  for (int i = 0 ; i < numberIntegers ; i++) {
781    int iColumn = integerVariable[i] ;
782    double djValue = direction*reducedCost[iColumn] ;
783    if (upper[iColumn]-lower[iColumn] > integerTolerance) {
784      if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap) {
785#ifdef COIN_HAS_CLP
786        // may just have been fixed before
787        if (clpSimplex) {
788          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
789#ifdef COIN_DEVELOP
790            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
791                   iColumn,clpSimplex->getColumnStatus(iColumn),
792                   djValue,gap,lower[iColumn],upper[iColumn]);
793#endif
794          } else {         
795            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
796                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
797          }
798        }
799#endif
800        solver->setColUpper(iColumn,lower[iColumn]) ;
801        numberFixed++ ;
802      } else if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap) {
803#ifdef COIN_HAS_CLP
804        // may just have been fixed before
805        if (clpSimplex) {
806          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
807#ifdef COIN_DEVELOP
808            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
809                   iColumn,clpSimplex->getColumnStatus(iColumn),
810                   djValue,gap,lower[iColumn],upper[iColumn]);
811#endif
812          } else {         
813            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
814                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
815          }
816        }
817#endif
818        solver->setColLower(iColumn,upper[iColumn]) ;
819        numberFixed++ ;
820      }
[916]821    }
822  }
[944]823  return numberFixed; 
[916]824}
Note: See TracBrowser for help on using the repository browser.