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

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

Corrected termination criterion.

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