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

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

more iterations at root

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