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

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

some changes e.g. DINS added and a bit of printing

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