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

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

compiler warnings, deterministic parallel and stability

File size: 25.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  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 = static_cast<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#ifdef DIVE_FIX_BINARY_VARIABLES
271  const double* reducedCost = solver->getReducedCost();
272#endif
273
274  int iteration = 0;
275  while(numberFractionalVariables) {
276    iteration++;
277
278    // select a fractional variable to bound
279    int bestColumn = -1;
280    int bestRound; // -1 rounds down, +1 rounds up
281    bool canRound = selectVariableToBranch(solver, newSolution, 
282                                           bestColumn, bestRound);
283
284    // if the solution is not trivially roundable, we don't try to round;
285    // if the solution is trivially roundable, we try to round. However,
286    // if the rounded solution is worse than the current incumbent,
287    // then we don't round and proceed normally. In this case, the
288    // bestColumn will be a trivially roundable variable
289    if(canRound) {
290      // check if by rounding all fractional variables
291      // we get a solution with an objective value
292      // better than the current best integer solution
293      double delta = 0.0;
294      for (int i=0; i<numberIntegers; i++) {
295        int iColumn = integerVariable[i];
296        double value=newSolution[iColumn];
297        if (fabs(floor(value+0.5)-value)>integerTolerance) {
298          assert(downLocks_[i]==0 || upLocks_[i]==0);
299          double obj = objective[iColumn];
300          if(downLocks_[i]==0 && upLocks_[i]==0) {
301            if(direction * obj >= 0.0)
302              delta += (floor(value) - value) * obj;
303            else
304              delta += (ceil(value) - value) * obj;
305          }
306          else if(downLocks_[i]==0)
307            delta += (floor(value) - value) * obj;
308          else
309            delta += (ceil(value) - value) * obj;
310        }
311      }
312      if(direction*(solver->getObjValue()+delta) < solutionValue) {
313#ifdef DIVE_DEBUG
314        nRoundFeasible++;
315#endif
316        // Round all the fractional variables
317        for (int i=0; i<numberIntegers; i++) {
318          int iColumn = integerVariable[i];
319          double value=newSolution[iColumn];
320          if (fabs(floor(value+0.5)-value)>integerTolerance) {
321            assert(downLocks_[i]==0 || upLocks_[i]==0);
322            if(downLocks_[i]==0 && upLocks_[i]==0) {
323              if(direction * objective[iColumn] >= 0.0)
324                newSolution[iColumn] = floor(value);
325              else
326                newSolution[iColumn] = ceil(value);
327            }
328            else if(downLocks_[i]==0)
329              newSolution[iColumn] = floor(value);
330            else
331              newSolution[iColumn] = ceil(value);
332          }
333        }
334        break;
335      }
336#ifdef DIVE_DEBUG
337      else
338        nRoundInfeasible++;
339#endif
340    }
341
342    // do reduced cost fixing
343#ifdef DIVE_DEBUG
344    int numberFixed = reducedCostFix(solver);
345    std::cout<<"numberReducedCostFixed = "<<numberFixed<<std::endl;
346#else
347    reducedCostFix(solver);
348#endif
349     
350    int numberAtBoundFixed = 0;
351#ifdef DIVE_FIX_BINARY_VARIABLES
352    // fix binary variables based on pseudo reduced cost
353    if(binVarIndex_.size()) {
354      int cnt = 0;
355      for (int j=0; j<(int)binVarIndex_.size(); j++) {
356        int iColumn1 = binVarIndex_[j];
357        double value = newSolution[iColumn1];
358        double maxPseudoReducedCost = 0.0;
359        if(fabs(value)<=integerTolerance &&
360           lower[iColumn1] != upper[iColumn1]) {
361#ifdef DIVE_DEBUG
362          std::cout<<"iColumn1 = "<<iColumn1<<", value = "<<value<<std::endl;
363#endif
364          int iRow = vbRowIndex_[j];
365          for (int k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {
366            int iColumn2 = column[k];
367#ifdef DIVE_DEBUG
368            std::cout<<"iColumn2 = "<<iColumn2<<std::endl;
369#endif
370            if(iColumn1 != iColumn2) {
371              double pseudoReducedCost = fabs(reducedCost[iColumn2] *
372                                              elementByRow[iColumn2] / 
373                                              elementByRow[iColumn1]);
374#ifdef DIVE_DEBUG
375              std::cout<<"reducedCost["<<iColumn2<<"] = "
376                       <<reducedCost[iColumn2]
377                       <<", elementByRow["<<iColumn2<<"] = "<<elementByRow[iColumn2]
378                       <<", elementByRow["<<iColumn1<<"] = "<<elementByRow[iColumn1]
379                       <<", pseudoRedCost = "<<pseudoReducedCost
380                       <<std::endl;
381#endif
382              if(pseudoReducedCost > maxPseudoReducedCost)
383                maxPseudoReducedCost = pseudoReducedCost;
384            }
385          }
386#ifdef DIVE_DEBUG
387          std::cout<<", maxPseudoRedCost = "<<maxPseudoReducedCost<<std::endl;
388#endif
389          candidate[cnt].var = iColumn1;
390          candidate[cnt++].pseudoRedCost = maxPseudoReducedCost;
391        }
392      }
393#ifdef DIVE_DEBUG
394      std::cout<<"candidates for rounding = "<<cnt<<std::endl;
395#endif
396      std::sort(candidate, candidate+cnt, compareBinaryVars);
397      for (int i=0; i<cnt; i++) {
398        int iColumn = candidate[i].var;
399        if (numberAtBoundFixed < maxNumberAtBoundToFix) {
400          columnFixed[numberAtBoundFixed] = iColumn;
401          originalBound[numberAtBoundFixed] = upper[iColumn];
402          fixedAtLowerBound[numberAtBoundFixed] = true;
403          solver->setColUpper(iColumn, lower[iColumn]);
404          numberAtBoundFixed++;
405          if(numberAtBoundFixed == maxNumberAtBoundToFix)
406            break;
407        }
408      }
409    }
410#endif
411
412    // fix other integer variables that are at their bounds
413    for (int i=0; i<numberIntegers; i++) {
414      int iColumn = integerVariable[i];
415      double value=newSolution[iColumn];
416      if(fabs(floor(value+0.5)-value)<=integerTolerance && 
417         numberAtBoundFixed < maxNumberAtBoundToFix) {
418        // fix the variable at one of its bounds
419        if (fabs(lower[iColumn]-value)<=integerTolerance &&
420            lower[iColumn] != upper[iColumn]) {
421          columnFixed[numberAtBoundFixed] = iColumn;
422          originalBound[numberAtBoundFixed] = upper[iColumn];
423          fixedAtLowerBound[numberAtBoundFixed] = true;
424          solver->setColUpper(iColumn, lower[iColumn]);
425          numberAtBoundFixed++;
426        }
427        else if(fabs(upper[iColumn]-value)<=integerTolerance &&
428            lower[iColumn] != upper[iColumn]) {
429          columnFixed[numberAtBoundFixed] = iColumn;
430          originalBound[numberAtBoundFixed] = lower[iColumn];
431          fixedAtLowerBound[numberAtBoundFixed] = false;
432          solver->setColLower(iColumn, upper[iColumn]);
433          numberAtBoundFixed++;
434        }
435        if(numberAtBoundFixed == maxNumberAtBoundToFix)
436          break;
437      }
438    }
439#ifdef DIVE_DEBUG
440    std::cout<<"numberAtBoundFixed = "<<numberAtBoundFixed<<std::endl;
441#endif
442
443    double originalBoundBestColumn;
444    if(bestColumn >= 0) {
445      if(bestRound < 0) {
446        originalBoundBestColumn = upper[bestColumn];
447        solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
448      }
449      else {
450        originalBoundBestColumn = lower[bestColumn];
451        solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
452      }
453    } else {
454      break;
455    }
456    int originalBestRound = bestRound;
457    int saveModelOptions = model_->specialOptions();
458    while (1) {
459
460      model_->setSpecialOptions(saveModelOptions|2048);
461      solver->resolve();
462      model_->setSpecialOptions(saveModelOptions);
463
464      if(!solver->isProvenOptimal()) {
465        if(numberAtBoundFixed > 0) {
466          // Remove the bound fix for variables that were at bounds
467          for(int i=0; i<numberAtBoundFixed; i++) {
468            int iColFixed = columnFixed[i];
469            if(fixedAtLowerBound[i])
470              solver->setColUpper(iColFixed, originalBound[i]);
471            else
472              solver->setColLower(iColFixed, originalBound[i]);
473          }
474          numberAtBoundFixed = 0;
475        }
476        else if(bestRound == originalBestRound) {
477          bestRound *= (-1);
478          if(bestRound < 0) {
479            solver->setColLower(bestColumn, originalBoundBestColumn);
480            solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
481          }
482          else {
483            solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
484            solver->setColUpper(bestColumn, originalBoundBestColumn);
485          }
486        }
487        else
488          break;
489      }
490      else
491        break;
492    }
493
494    if(!solver->isProvenOptimal() || 
495       direction*solver->getObjValue() >= solutionValue) {
496#ifdef DIVE_DEBUG
497      reasonToStop = 1;
498#endif
499      break;
500    }
501
502    if(iteration > maxIterations_) {
503#ifdef DIVE_DEBUG
504      reasonToStop = 2;
505#endif
506      break;
507    }
508
509    if(CoinCpuTime()-time1 > maxTime_) {
510#ifdef DIVE_DEBUG
511      reasonToStop = 3;
512#endif
513      break;
514    }
515
516    numberSimplexIterations+=solver->getIterationCount();
517    if(numberSimplexIterations > maxSimplexIterations) {
518#ifdef DIVE_DEBUG
519      reasonToStop = 4;
520#endif
521      // also switch off
522#ifdef CLP_INVESTIGATE
523      printf("switching off diving as too many iterations %d, %d allowed\n",
524             numberSimplexIterations,maxSimplexIterations);
525#endif
526      when_=0;
527      break;
528    }
529
530    if (solver->getIterationCount()>1000&&iteration>3) {
531#ifdef DIVE_DEBUG
532      reasonToStop = 5;
533#endif
534      // also switch off
535#ifdef CLP_INVESTIGATE
536      printf("switching off diving one iteration took %d iterations (total %d)\n",
537             solver->getIterationCount(),numberSimplexIterations);
538#endif
539      when_=0;
540      break;
541    }
542
543    memcpy(newSolution,solution,numberColumns*sizeof(double));
544    numberFractionalVariables = 0;
545    for (int i=0; i<numberIntegers; i++) {
546      int iColumn = integerVariable[i];
547      double value=newSolution[iColumn];
548      if (fabs(floor(value+0.5)-value)>integerTolerance) {
549        numberFractionalVariables++;
550      }
551    }
552
553  }
554
555
556  double * rowActivity = new double[numberRows];
557  memset(rowActivity,0,numberRows*sizeof(double));
558
559  // re-compute new solution value
560  double objOffset=0.0;
561  solver->getDblParam(OsiObjOffset,objOffset);
562  newSolutionValue = -objOffset;
563  for (int i=0 ; i<numberColumns ; i++ )
564    newSolutionValue += objective[i]*newSolution[i];
565  newSolutionValue *= direction;
566    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
567  if (newSolutionValue<solutionValue) {
568    // paranoid check
569    memset(rowActivity,0,numberRows*sizeof(double));
570    for (int i=0;i<numberColumns;i++) {
571      int j;
572      double value = newSolution[i];
573      if (value) {
574        for (j=columnStart[i];
575             j<columnStart[i]+columnLength[i];j++) {
576          int iRow=row[j];
577          rowActivity[iRow] += value*element[j];
578        }
579      }
580    }
581    // check was approximately feasible
582    bool feasible=true;
583    for (int i=0;i<numberRows;i++) {
584      if(rowActivity[i]<rowLower[i]) {
585        if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
586          feasible = false;
587      } else if(rowActivity[i]>rowUpper[i]) {
588        if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
589          feasible = false;
590      }
591    }
592    for (int i=0; i<numberIntegers; i++) {
593      int iColumn = integerVariable[i];
594      double value=newSolution[iColumn];
595      if (fabs(floor(value+0.5)-value)>integerTolerance) {
596        feasible = false;
597        break;
598      }
599    }
600    if (feasible) {
601      // new solution
602      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
603      solutionValue = newSolutionValue;
604      //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
605      returnCode=1;
606    } else {
607      // Can easily happen
608      //printf("Debug CbcHeuristicDive giving bad solution\n");
609    }
610  }
611
612#ifdef DIVE_DEBUG
613  std::cout<<"nRoundInfeasible = "<<nRoundInfeasible
614           <<", nRoundFeasible = "<<nRoundFeasible
615           <<", returnCode = "<<returnCode
616           <<", reasonToStop = "<<reasonToStop
617           <<", simplexIts = "<<numberSimplexIterations
618           <<", iterations = "<<iteration<<std::endl;
619#endif
620
621  delete [] newSolution;
622  delete [] columnFixed;
623  delete [] originalBound;
624  delete [] fixedAtLowerBound;
625  delete [] candidate;
626  delete [] rowActivity;
627  delete solver;
628  return returnCode;
629}
630
631// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
632void 
633CbcHeuristicDive::validate() 
634{
635  if (model_&&when()<10) {
636    if (model_->numberIntegers()!=
637        model_->numberObjects()&&(model_->numberObjects()||
638                                  (model_->specialOptions()&1024)==0))
639      setWhen(0);
640  }
641
642  int numberIntegers = model_->numberIntegers();
643  const int * integerVariable = model_->integerVariable();
644  delete [] downLocks_;
645  delete [] upLocks_;
646  downLocks_ = new unsigned short [numberIntegers];
647  upLocks_ = new unsigned short [numberIntegers];
648  // Column copy
649  const double * element = matrix_.getElements();
650  const int * row = matrix_.getIndices();
651  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
652  const int * columnLength = matrix_.getVectorLengths();
653  const double * rowLower = model_->solver()->getRowLower();
654  const double * rowUpper = model_->solver()->getRowUpper();
655  for (int i=0;i<numberIntegers;i++) {
656    int iColumn = integerVariable[i];
657    int down=0;
658    int up=0;
659    if (columnLength[iColumn]>65535) {
660      setWhen(0);
661      break; // unlikely to work
662    }
663    for (CoinBigIndex j=columnStart[iColumn];
664         j<columnStart[iColumn]+columnLength[iColumn];j++) {
665      int iRow=row[j];
666      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
667        up++;
668        down++;
669      } else if (element[j]>0.0) {
670        if (rowUpper[iRow]<1.0e20)
671          up++;
672        else
673          down++;
674      } else {
675        if (rowLower[iRow]>-1.0e20)
676          up++;
677        else
678          down++;
679      }
680    }
681    downLocks_[i] = static_cast<unsigned short> (down);
682    upLocks_[i] = static_cast<unsigned short> (up);
683  }
684
685#ifdef DIVE_FIX_BINARY_VARIABLES
686  selectBinaryVariables();
687#endif
688}
689
690// Select candidate binary variables for fixing
691void
692CbcHeuristicDive::selectBinaryVariables()
693{
694  // Row copy
695  const double * elementByRow = matrixByRow_.getElements();
696  const int * column = matrixByRow_.getIndices();
697  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
698  const int * rowLength = matrixByRow_.getVectorLengths();
699
700  const int numberRows = matrixByRow_.getNumRows();
701  const int numberCols = matrixByRow_.getNumCols();
702
703  const double * lower = model_->solver()->getColLower();
704  const double * upper = model_->solver()->getColUpper();
705  const double * rowLower = model_->solver()->getRowLower();
706  const double * rowUpper = model_->solver()->getRowUpper();
707
708  //  const char * integerType = model_->integerType();
709 
710
711  //  const int numberIntegers = model_->numberIntegers();
712  //  const int * integerVariable = model_->integerVariable();
713  const double * objective = model_->solver()->getObjCoefficients();
714
715  // vector to store the row number of variable bound rows
716  int* rowIndexes = new int [numberCols];
717  memset(rowIndexes, -1, numberCols*sizeof(int));
718
719  for(int i=0; i<numberRows; i++) {
720    int positiveBinary = -1;
721    int negativeBinary = -1;
722    int nPositiveOther = 0;
723    int nNegativeOther = 0;
724    for (int k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
725      int iColumn = column[k];
726      if(model_->solver()->isInteger(iColumn) &&
727         lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
728         objective[iColumn] == 0.0 &&
729         elementByRow[k] > 0.0 &&
730         positiveBinary < 0)
731        positiveBinary = iColumn;
732      else if(model_->solver()->isInteger(iColumn) &&
733              lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
734              objective[iColumn] == 0.0 &&
735              elementByRow[k] < 0.0 &&
736              negativeBinary < 0)
737        negativeBinary = iColumn;
738      else if((elementByRow[k]> 0.0 &&
739               lower[iColumn] >= 0.0) ||
740              (elementByRow[k] < 0.0 &&
741               upper[iColumn] <= 0.0))
742        nPositiveOther++;
743      else if((elementByRow[k]> 0.0 &&
744               lower[iColumn] <= 0.0) ||
745              (elementByRow[k] < 0.0 &&
746               upper[iColumn] >= 0.0))
747        nNegativeOther++;
748      if(nPositiveOther > 0 && nNegativeOther > 0)
749        break;
750    }
751    int binVar = -1;
752    if(positiveBinary >= 0 &&
753       (negativeBinary >= 0 || nNegativeOther > 0) &&
754       nPositiveOther == 0 &&
755       rowLower[i] == 0.0 &&
756       rowUpper[i] > 0.0)
757      binVar = positiveBinary;
758    else if(negativeBinary >= 0 &&
759            (positiveBinary >= 0 || nPositiveOther > 0) &&
760            nNegativeOther == 0 &&
761            rowLower[i] < 0.0 &&
762            rowUpper[i] == 0.0)
763      binVar = negativeBinary;
764    if(binVar >= 0) {
765      if(rowIndexes[binVar] == -1)
766        rowIndexes[binVar] = i;
767      else if(rowIndexes[binVar] >= 0)
768        rowIndexes[binVar] = -2;
769    }
770  }
771
772  for(int j=0; j<numberCols; j++) {
773    if(rowIndexes[j] >= 0) {
774      binVarIndex_.push_back(j);
775      vbRowIndex_.push_back(rowIndexes[j]);
776    }
777  }
778
779#ifdef DIVE_DEBUG
780  std::cout<<"number vub Binary = "<<binVarIndex_.size()<<std::endl;
781#endif
782
783  delete [] rowIndexes;
784   
785}
786
787/*
788  Perform reduced cost fixing on integer variables.
789
790  The variables in question are already nonbasic at bound. We're just nailing
791  down the current situation.
792*/
793
794int CbcHeuristicDive::reducedCostFix (OsiSolverInterface* solver)
795
796{
797  return 0; // temp
798#if 0
799  if(!solverCharacteristics_->reducedCostsAccurate())
800    return 0; //NLP
801#endif
802  double cutoff = model_->getCutoff() ;
803#ifdef DIVE_DEBUG
804  std::cout<<"cutoff = "<<cutoff<<std::endl;
805#endif
806  double direction = solver->getObjSense() ;
807  double gap = cutoff - solver->getObjValue()*direction ;
808  double tolerance;
809  solver->getDblParam(OsiDualTolerance,tolerance) ;
810  if (gap<=0.0)
811    gap = tolerance; //return 0;
812  gap += 100.0*tolerance;
813  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
814
815  const double *lower = solver->getColLower() ;
816  const double *upper = solver->getColUpper() ;
817  const double *solution = solver->getColSolution() ;
818  const double *reducedCost = solver->getReducedCost() ;
819
820  int numberIntegers = model_->numberIntegers();
821  const int * integerVariable = model_->integerVariable();
822
823  int numberFixed = 0 ;
824
825# ifdef COIN_HAS_CLP
826  OsiClpSolverInterface * clpSolver
827    = dynamic_cast<OsiClpSolverInterface *> (solver);
828  ClpSimplex * clpSimplex=NULL;
829  if (clpSolver) 
830    clpSimplex = clpSolver->getModelPtr();
831# endif
832  for (int i = 0 ; i < numberIntegers ; i++) {
833    int iColumn = integerVariable[i] ;
834    double djValue = direction*reducedCost[iColumn] ;
835    if (upper[iColumn]-lower[iColumn] > integerTolerance) {
836      if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap) {
837#ifdef COIN_HAS_CLP
838        // may just have been fixed before
839        if (clpSimplex) {
840          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
841#ifdef COIN_DEVELOP
842            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
843                   iColumn,clpSimplex->getColumnStatus(iColumn),
844                   djValue,gap,lower[iColumn],upper[iColumn]);
845#endif
846          } else {         
847            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
848                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
849          }
850        }
851#endif
852        solver->setColUpper(iColumn,lower[iColumn]) ;
853        numberFixed++ ;
854      } else if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap) {
855#ifdef COIN_HAS_CLP
856        // may just have been fixed before
857        if (clpSimplex) {
858          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
859#ifdef COIN_DEVELOP
860            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
861                   iColumn,clpSimplex->getColumnStatus(iColumn),
862                   djValue,gap,lower[iColumn],upper[iColumn]);
863#endif
864          } else {         
865            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
866                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
867          }
868        }
869#endif
870        solver->setColLower(iColumn,upper[iColumn]) ;
871        numberFixed++ ;
872      }
873    }
874  }
875  return numberFixed; 
876}
Note: See TracBrowser for help on using the repository browser.