source: stable/2.4/Cbc/src/CbcHeuristicDive.cpp @ 1271

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

Creating new stable branch 2.4 from trunk (rev 1270)

File size: 31.6 KB
Line 
1/* $Id: CbcHeuristicDive.cpp 1240 2009-10-02 18:41:44Z forrest $ */
2// Copyright (C) 2008, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8
9#include "CbcHeuristicDive.hpp"
10#include "CbcStrategy.hpp"
11#include "OsiAuxInfo.hpp"
12#include  "CoinTime.hpp"
13
14#ifdef COIN_HAS_CLP
15#include "OsiClpSolverInterface.hpp"
16#endif
17
18//#define DIVE_FIX_BINARY_VARIABLES
19//#define DIVE_DEBUG
20
21// Default Constructor
22CbcHeuristicDive::CbcHeuristicDive() 
23  :CbcHeuristic()
24{
25  // matrix and row copy will automatically be empty
26  downLocks_ =NULL;
27  upLocks_ = NULL;
28  downArray_ =NULL;
29  upArray_ = NULL;
30  percentageToFix_ = 0.2;
31  maxIterations_ = 100;
32  maxSimplexIterations_ = 10000;
33  maxSimplexIterationsAtRoot_ = 1000000;
34  maxTime_ = 600;
35  whereFrom_=255-2-16+256;
36}
37
38// Constructor from model
39CbcHeuristicDive::CbcHeuristicDive(CbcModel & model)
40  :CbcHeuristic(model)
41{
42  downLocks_ =NULL;
43  upLocks_ = NULL;
44  downArray_ =NULL;
45  upArray_ = NULL;
46  // Get a copy of original matrix
47  assert(model.solver());
48  // model may have empty matrix - wait until setModel
49  const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
50  if (matrix) {
51    matrix_ = *matrix;
52    matrixByRow_ = *model.solver()->getMatrixByRow();
53    validate();
54  }
55  percentageToFix_ = 0.2;
56  maxIterations_ = 100;
57  maxSimplexIterations_ = 10000;
58  maxSimplexIterationsAtRoot_ = 1000000;
59  maxTime_ = 600;
60  whereFrom_=255-2-16+256;
61}
62
63// Destructor
64CbcHeuristicDive::~CbcHeuristicDive ()
65{
66  delete [] downLocks_;
67  delete [] upLocks_;
68  assert (!downArray_);
69}
70
71// Create C++ lines to get to current state
72void 
73CbcHeuristicDive::generateCpp( FILE * fp, const char * heuristic) 
74{
75  // hard coded as CbcHeuristic virtual
76  CbcHeuristic::generateCpp(fp,heuristic);
77  if (percentageToFix_!=0.2)
78    fprintf(fp,"3  %s.setPercentageToFix(%.f);\n",heuristic,percentageToFix_);
79  else
80    fprintf(fp,"4  %s.setPercentageToFix(%.f);\n",heuristic,percentageToFix_);
81  if (maxIterations_!=100)
82    fprintf(fp,"3  %s.setMaxIterations(%d);\n",heuristic,maxIterations_);
83  else
84    fprintf(fp,"4  %s.setMaxIterations(%d);\n",heuristic,maxIterations_);
85  if (maxSimplexIterations_!=10000)
86    fprintf(fp,"3  %s.setMaxSimplexIterations(%d);\n",heuristic,maxSimplexIterations_);
87  else
88    fprintf(fp,"4  %s.setMaxSimplexIterations(%d);\n",heuristic,maxSimplexIterations_);
89  if (maxTime_!=600)
90    fprintf(fp,"3  %s.setMaxTime(%.2f);\n",heuristic,maxTime_);
91  else
92    fprintf(fp,"4  %s.setMaxTime(%.2f);\n",heuristic,maxTime_);
93}
94
95// Copy constructor
96CbcHeuristicDive::CbcHeuristicDive(const CbcHeuristicDive & rhs)
97:
98  CbcHeuristic(rhs),
99  matrix_(rhs.matrix_),
100  matrixByRow_(rhs.matrixByRow_),
101  percentageToFix_(rhs.percentageToFix_),
102  maxIterations_(rhs.maxIterations_),
103  maxSimplexIterations_(rhs.maxSimplexIterations_),
104  maxSimplexIterationsAtRoot_(rhs.maxSimplexIterationsAtRoot_),
105  maxTime_(rhs.maxTime_)
106{
107  downArray_ =NULL;
108  upArray_ = NULL;
109  if (rhs.downLocks_) {
110    int numberIntegers = model_->numberIntegers();
111    downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
112    upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
113  } else {
114    downLocks_ = NULL;
115    upLocks_ = NULL;
116  }
117}
118
119// Assignment operator
120CbcHeuristicDive & 
121CbcHeuristicDive::operator=( const CbcHeuristicDive& rhs)
122{
123  if (this!=&rhs) {
124    CbcHeuristic::operator=(rhs);
125    matrix_ = rhs.matrix_;
126    matrixByRow_ = rhs.matrixByRow_;
127    percentageToFix_ = rhs.percentageToFix_;
128    maxIterations_ = rhs.maxIterations_;
129    maxSimplexIterations_ = rhs.maxSimplexIterations_;
130    maxSimplexIterationsAtRoot_ = rhs.maxSimplexIterationsAtRoot_;
131    maxTime_ = rhs.maxTime_;
132    delete [] downLocks_;
133    delete [] upLocks_;
134    if (rhs.downLocks_) {
135      int numberIntegers = model_->numberIntegers();
136      downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
137      upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
138    } else {
139      downLocks_ = NULL;
140      upLocks_ = NULL;
141    }
142  }
143  return *this;
144}
145
146// Resets stuff if model changes
147void 
148CbcHeuristicDive::resetModel(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  // model may have empty matrix - wait until setModel
155  if (matrix) {
156    matrix_ = *matrix;
157    matrixByRow_ = *model->solver()->getMatrixByRow();
158    validate();
159  }
160}
161
162// update model
163void CbcHeuristicDive::setModel(CbcModel * model)
164{
165  model_ = model;
166  assert(model_->solver());
167  // Get a copy of original matrix
168  const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
169  if (matrix) {
170    matrix_ = *matrix;
171    matrixByRow_ = *model->solver()->getMatrixByRow();
172    // make sure model okay for heuristic
173    validate();
174  }
175}
176
177bool CbcHeuristicDive::canHeuristicRun()
178{
179  return shouldHeurRun_randomChoice();
180}
181
182inline bool compareBinaryVars(const PseudoReducedCost obj1,
183                              const PseudoReducedCost obj2)
184{
185  return obj1.pseudoRedCost > obj2.pseudoRedCost;
186}
187
188// See if diving will give better solution
189// Sets value of solution
190// Returns 1 if solution, 0 if not
191int
192CbcHeuristicDive::solution(double & solutionValue,
193                           double * betterSolution)
194{
195#ifdef DIVE_DEBUG
196  std::cout<<"solutionValue = "<<solutionValue<<std::endl;
197#endif
198  ++numCouldRun_;
199
200  // test if the heuristic can run
201  if(!canHeuristicRun())
202    return 0;
203
204#if 0
205  // See if to do
206  if (!when()||(when()%10==1&&model_->phase()!=1)||
207      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
208    return 0; // switched off
209#endif
210
211#ifdef DIVE_DEBUG
212  int nRoundInfeasible = 0;
213  int nRoundFeasible = 0;
214  int reasonToStop = 0;
215#endif
216
217  double time1 = CoinCpuTime();
218  int numberSimplexIterations=0;
219  int maxSimplexIterations= (model_->getNodeCount()) ? maxSimplexIterations_
220    : maxSimplexIterationsAtRoot_;
221
222  OsiSolverInterface * solver = cloneBut(6); // was model_->solver()->clone();
223# ifdef COIN_HAS_CLP
224  OsiClpSolverInterface * clpSolver
225    = dynamic_cast<OsiClpSolverInterface *> (solver);
226  if (clpSolver) {
227    // say give up easily
228    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
229    clpSimplex->setMoreSpecialOptions(clpSimplex->moreSpecialOptions()|64);
230  }
231# endif
232  const double * lower = solver->getColLower();
233  const double * upper = solver->getColUpper();
234  const double * rowLower = solver->getRowLower();
235  const double * rowUpper = solver->getRowUpper();
236  const double * solution = solver->getColSolution();
237  const double * objective = solver->getObjCoefficients();
238  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
239  double primalTolerance;
240  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
241
242  int numberRows = matrix_.getNumRows();
243  assert (numberRows<=solver->getNumRows());
244  int numberIntegers = model_->numberIntegers();
245  const int * integerVariable = model_->integerVariable();
246  double direction = solver->getObjSense(); // 1 for min, -1 for max
247  double newSolutionValue = direction*solver->getObjValue();
248  int returnCode = 0;
249  // Column copy
250  const double * element = matrix_.getElements();
251  const int * row = matrix_.getIndices();
252  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
253  const int * columnLength = matrix_.getVectorLengths();
254#ifdef DIVE_FIX_BINARY_VARIABLES
255  // Row copy
256  const double * elementByRow = matrixByRow_.getElements();
257  const int * column = matrixByRow_.getIndices();
258  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
259  const int * rowLength = matrixByRow_.getVectorLengths();
260#endif
261
262  // Get solution array for heuristic solution
263  int numberColumns = solver->getNumCols();
264  double * newSolution = new double [numberColumns];
265  memcpy(newSolution,solution,numberColumns*sizeof(double));
266
267  // vectors to store the latest variables fixed at their bounds
268  int* columnFixed = new int [numberIntegers];
269  double* originalBound = new double [numberIntegers];
270  bool * fixedAtLowerBound = new bool [numberIntegers];
271  PseudoReducedCost * candidate = new PseudoReducedCost [numberIntegers];
272  double * random = new double [numberIntegers];
273
274  int maxNumberAtBoundToFix = static_cast<int> (floor(percentageToFix_ * numberIntegers));
275
276  // count how many fractional variables
277  int numberFractionalVariables = 0;
278  for (int i=0; i<numberIntegers; i++) {
279    random[i]=randomNumberGenerator_.randomDouble()+0.3;
280    int iColumn = integerVariable[i];
281    double value=newSolution[iColumn];
282    if (fabs(floor(value+0.5)-value)>integerTolerance) {
283      numberFractionalVariables++;
284    }
285  }
286
287  const double* reducedCost = NULL;
288  // See if not NLP
289  if(model_->solverCharacteristics()->reducedCostsAccurate())
290    reducedCost = solver->getReducedCost();
291
292  int iteration = 0;
293  while(numberFractionalVariables) {
294    iteration++;
295   
296    // initialize any data
297    initializeData();
298
299    // select a fractional variable to bound
300    int bestColumn = -1;
301    int bestRound; // -1 rounds down, +1 rounds up
302    bool canRound = selectVariableToBranch(solver, newSolution, 
303                                           bestColumn, bestRound);
304
305    // if the solution is not trivially roundable, we don't try to round;
306    // if the solution is trivially roundable, we try to round. However,
307    // if the rounded solution is worse than the current incumbent,
308    // then we don't round and proceed normally. In this case, the
309    // bestColumn will be a trivially roundable variable
310    if(canRound) {
311      // check if by rounding all fractional variables
312      // we get a solution with an objective value
313      // better than the current best integer solution
314      double delta = 0.0;
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          double obj = objective[iColumn];
321          if(downLocks_[i]==0 && upLocks_[i]==0) {
322            if(direction * obj >= 0.0)
323              delta += (floor(value) - value) * obj;
324            else
325              delta += (ceil(value) - value) * obj;
326          }
327          else if(downLocks_[i]==0)
328            delta += (floor(value) - value) * obj;
329          else
330            delta += (ceil(value) - value) * obj;
331        }
332      }
333      if(direction*(solver->getObjValue()+delta) < solutionValue) {
334#ifdef DIVE_DEBUG
335        nRoundFeasible++;
336#endif
337        // Round all the fractional variables
338        for (int i=0; i<numberIntegers; i++) {
339          int iColumn = integerVariable[i];
340          double value=newSolution[iColumn];
341          if (fabs(floor(value+0.5)-value)>integerTolerance) {
342            assert(downLocks_[i]==0 || upLocks_[i]==0);
343            if(downLocks_[i]==0 && upLocks_[i]==0) {
344              if(direction * objective[iColumn] >= 0.0)
345                newSolution[iColumn] = floor(value);
346              else
347                newSolution[iColumn] = ceil(value);
348            }
349            else if(downLocks_[i]==0)
350              newSolution[iColumn] = floor(value);
351            else
352              newSolution[iColumn] = ceil(value);
353          }
354        }
355        break;
356      }
357#ifdef DIVE_DEBUG
358      else
359        nRoundInfeasible++;
360#endif
361    }
362
363    // do reduced cost fixing
364#ifdef DIVE_DEBUG
365    int numberFixed = reducedCostFix(solver);
366    std::cout<<"numberReducedCostFixed = "<<numberFixed<<std::endl;
367#else
368    reducedCostFix(solver);
369#endif
370     
371    int numberAtBoundFixed = 0;
372#ifdef DIVE_FIX_BINARY_VARIABLES
373    // fix binary variables based on pseudo reduced cost
374    if(binVarIndex_.size()) {
375      int cnt = 0;
376      int n = static_cast<int>(binVarIndex_.size());
377      for (int j=0; j<n; j++) {
378        int iColumn1 = binVarIndex_[j];
379        double value = newSolution[iColumn1];
380        if(fabs(value)<=integerTolerance &&
381           lower[iColumn1] != upper[iColumn1]) {
382          double maxPseudoReducedCost = 0.0;
383#ifdef DIVE_DEBUG
384          std::cout<<"iColumn1 = "<<iColumn1<<", value = "<<value<<std::endl;
385#endif
386          int iRow = vbRowIndex_[j];
387          double chosenValue=0.0;
388          for (int k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {
389            int iColumn2 = column[k];
390#ifdef DIVE_DEBUG
391            std::cout<<"iColumn2 = "<<iColumn2<<std::endl;
392#endif
393            if(iColumn1 != iColumn2) {
394              double pseudoReducedCost = fabs(reducedCost[iColumn2] *
395                                              elementByRow[k]);
396#ifdef DIVE_DEBUG
397              int k2;
398              for (k2=rowStart[iRow];k2<rowStart[iRow]+rowLength[iRow];k2++) {
399                if (column[k2]==iColumn1)
400                  break;
401              }
402              std::cout<<"reducedCost["<<iColumn2<<"] = "
403                       <<reducedCost[iColumn2]
404                       <<", elementByRow["<<iColumn2<<"] = "<<elementByRow[k]
405                       <<", elementByRow["<<iColumn1<<"] = "<<elementByRow[k2]
406                       <<", pseudoRedCost = "<<pseudoReducedCost
407                       <<std::endl;
408#endif
409              if(pseudoReducedCost > maxPseudoReducedCost)
410                maxPseudoReducedCost = pseudoReducedCost;
411            } else {
412              // save value
413              chosenValue = fabs(elementByRow[k]);
414            }
415          }
416          assert (chosenValue);
417          maxPseudoReducedCost /= chosenValue;
418#ifdef DIVE_DEBUG
419          std::cout<<", maxPseudoRedCost = "<<maxPseudoReducedCost<<std::endl;
420#endif
421          candidate[cnt].var = iColumn1;
422          candidate[cnt++].pseudoRedCost = maxPseudoReducedCost;
423        }
424      }
425#ifdef DIVE_DEBUG
426      std::cout<<"candidates for rounding = "<<cnt<<std::endl;
427#endif
428      std::sort(candidate, candidate+cnt, compareBinaryVars);
429      for (int i=0; i<cnt; i++) {
430        int iColumn = candidate[i].var;
431        if (numberAtBoundFixed < maxNumberAtBoundToFix) {
432          columnFixed[numberAtBoundFixed] = iColumn;
433          originalBound[numberAtBoundFixed] = upper[iColumn];
434          fixedAtLowerBound[numberAtBoundFixed] = true;
435          solver->setColUpper(iColumn, lower[iColumn]);
436          numberAtBoundFixed++;
437          if(numberAtBoundFixed == maxNumberAtBoundToFix)
438            break;
439        }
440      }
441    }
442#endif
443
444    // fix other integer variables that are at their bounds
445    int cnt=0;
446#ifdef GAP
447    double gap=1.0e30;
448#endif
449    if (reducedCost&&true) {
450#if 1
451      cnt=fixOtherVariables(solver,solution,candidate,random);
452#else
453#ifdef GAP
454      double cutoff = model_->getCutoff() ; 
455      if (cutoff<1.0e20&&false) {
456        double direction = solver->getObjSense() ;
457        gap = cutoff - solver->getObjValue()*direction ;
458        gap *= 0.1; // Fix more if plausible
459        double tolerance;
460        solver->getDblParam(OsiDualTolerance,tolerance) ;
461        if (gap<=0.0)
462          gap = tolerance;
463        gap += 100.0*tolerance;
464      }
465      int nOverGap=0;
466#endif
467      int numberFree=0;
468      int numberFixed=0;
469      for (int i=0; i<numberIntegers; i++) {
470        int iColumn = integerVariable[i];
471        if (upper[iColumn]>lower[iColumn]) {
472          numberFree++;
473          double value=newSolution[iColumn];
474          if(fabs(floor(value+0.5)-value)<=integerTolerance) {
475            candidate[cnt].var = iColumn;
476            candidate[cnt++].pseudoRedCost = 
477              fabs(reducedCost[iColumn]*random[i]);
478#ifdef GAP
479            if (fabs(reducedCost[iColumn])>gap)
480              nOverGap++;
481#endif
482          }
483        } else {
484          numberFixed++;
485        }
486      }
487#ifdef GAP
488      int nLeft = maxNumberAtBoundToFix-numberAtBoundFixed;
489#ifdef CLP_INVESTIGATE
490      printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
491             cutoff,solver->getObjValue(),nOverGap,numberFree,numberFixed);
492#endif
493      if (nOverGap>nLeft&&true) {
494        nOverGap = CoinMin(nOverGap,nLeft+maxNumberAtBoundToFix/2);
495        maxNumberAtBoundToFix += nOverGap-nLeft;
496      }
497#else
498#ifdef CLP_INVESTIGATE
499      printf("cutoff %g obj %g - %d free, %d fixed\n",
500             model_->getCutoff(),solver->getObjValue(),numberFree,numberFixed);
501#endif
502#endif
503#endif
504    } else {
505      for (int i=0; i<numberIntegers; i++) {
506        int iColumn = integerVariable[i];
507        if (upper[iColumn]>lower[iColumn]) {
508          double value=newSolution[iColumn];
509          if(fabs(floor(value+0.5)-value)<=integerTolerance) {
510            candidate[cnt].var = iColumn;
511            candidate[cnt++].pseudoRedCost = numberIntegers-i;
512          }
513        }
514      }
515    }
516    std::sort(candidate, candidate+cnt, compareBinaryVars);
517    for (int i=0; i<cnt; i++) {
518      int iColumn = candidate[i].var;
519      if (upper[iColumn]>lower[iColumn]) {
520        double value=newSolution[iColumn];
521        if(fabs(floor(value+0.5)-value)<=integerTolerance && 
522           numberAtBoundFixed < maxNumberAtBoundToFix) {
523          // fix the variable at one of its bounds
524          if (fabs(lower[iColumn]-value)<=integerTolerance) {
525            columnFixed[numberAtBoundFixed] = iColumn;
526            originalBound[numberAtBoundFixed] = upper[iColumn];
527            fixedAtLowerBound[numberAtBoundFixed] = true;
528            solver->setColUpper(iColumn, lower[iColumn]);
529            numberAtBoundFixed++;
530          }
531          else if(fabs(upper[iColumn]-value)<=integerTolerance) {
532            columnFixed[numberAtBoundFixed] = iColumn;
533            originalBound[numberAtBoundFixed] = lower[iColumn];
534            fixedAtLowerBound[numberAtBoundFixed] = false;
535            solver->setColLower(iColumn, upper[iColumn]);
536            numberAtBoundFixed++;
537          }
538          if(numberAtBoundFixed == maxNumberAtBoundToFix)
539            break;
540        }
541      }
542    }
543#ifdef DIVE_DEBUG
544    std::cout<<"numberAtBoundFixed = "<<numberAtBoundFixed<<std::endl;
545#endif
546
547    double originalBoundBestColumn;
548    if(bestColumn >= 0) {
549      if(bestRound < 0) {
550        originalBoundBestColumn = upper[bestColumn];
551        solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
552      }
553      else {
554        originalBoundBestColumn = lower[bestColumn];
555        solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
556      }
557    } else {
558      break;
559    }
560    int originalBestRound = bestRound;
561    int saveModelOptions = model_->specialOptions();
562    while (1) {
563
564      model_->setSpecialOptions(saveModelOptions|2048);
565      solver->resolve();
566      model_->setSpecialOptions(saveModelOptions);
567      if(!solver->isAbandoned()) {
568        numberSimplexIterations+=solver->getIterationCount();
569      } else {
570        numberSimplexIterations=maxSimplexIterations+1;
571        break;
572      }
573     
574      if(!solver->isProvenOptimal()) {
575        if(numberAtBoundFixed > 0) {
576          // Remove the bound fix for variables that were at bounds
577          for(int i=0; i<numberAtBoundFixed; i++) {
578            int iColFixed = columnFixed[i];
579            if(fixedAtLowerBound[i])
580              solver->setColUpper(iColFixed, originalBound[i]);
581            else
582              solver->setColLower(iColFixed, originalBound[i]);
583          }
584          numberAtBoundFixed = 0;
585        }
586        else if(bestRound == originalBestRound) {
587          bestRound *= (-1);
588          if(bestRound < 0) {
589            solver->setColLower(bestColumn, originalBoundBestColumn);
590            solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
591          }
592          else {
593            solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
594            solver->setColUpper(bestColumn, originalBoundBestColumn);
595          }
596        }
597        else
598          break;
599      }
600      else
601        break;
602    }
603
604    if(!solver->isProvenOptimal() || 
605       direction*solver->getObjValue() >= solutionValue) {
606#ifdef DIVE_DEBUG
607      reasonToStop = 1;
608#endif
609      break;
610    }
611
612    if(iteration > maxIterations_) {
613#ifdef DIVE_DEBUG
614      reasonToStop = 2;
615#endif
616      break;
617    }
618
619    if(CoinCpuTime()-time1 > maxTime_) {
620#ifdef DIVE_DEBUG
621      reasonToStop = 3;
622#endif
623      break;
624    }
625
626    if(numberSimplexIterations > maxSimplexIterations) {
627#ifdef DIVE_DEBUG
628      reasonToStop = 4;
629#endif
630      // also switch off
631#ifdef CLP_INVESTIGATE
632      printf("switching off diving as too many iterations %d, %d allowed\n",
633             numberSimplexIterations,maxSimplexIterations);
634#endif
635      when_=0;
636      break;
637    }
638
639    if (solver->getIterationCount()>1000&&iteration>3) {
640#ifdef DIVE_DEBUG
641      reasonToStop = 5;
642#endif
643      // also switch off
644#ifdef CLP_INVESTIGATE
645      printf("switching off diving one iteration took %d iterations (total %d)\n",
646             solver->getIterationCount(),numberSimplexIterations);
647#endif
648      when_=0;
649      break;
650    }
651
652    memcpy(newSolution,solution,numberColumns*sizeof(double));
653    numberFractionalVariables = 0;
654    for (int i=0; i<numberIntegers; i++) {
655      int iColumn = integerVariable[i];
656      double value=newSolution[iColumn];
657      if (fabs(floor(value+0.5)-value)>integerTolerance) {
658        numberFractionalVariables++;
659      }
660    }
661
662  }
663
664
665  double * rowActivity = new double[numberRows];
666  memset(rowActivity,0,numberRows*sizeof(double));
667
668  // re-compute new solution value
669  double objOffset=0.0;
670  solver->getDblParam(OsiObjOffset,objOffset);
671  newSolutionValue = -objOffset;
672  for (int i=0 ; i<numberColumns ; i++ )
673    newSolutionValue += objective[i]*newSolution[i];
674  newSolutionValue *= direction;
675    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
676  if (newSolutionValue<solutionValue) {
677    // paranoid check
678    memset(rowActivity,0,numberRows*sizeof(double));
679    for (int i=0;i<numberColumns;i++) {
680      int j;
681      double value = newSolution[i];
682      if (value) {
683        for (j=columnStart[i];
684             j<columnStart[i]+columnLength[i];j++) {
685          int iRow=row[j];
686          rowActivity[iRow] += value*element[j];
687        }
688      }
689    }
690    // check was approximately feasible
691    bool feasible=true;
692    for (int i=0;i<numberRows;i++) {
693      if(rowActivity[i]<rowLower[i]) {
694        if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
695          feasible = false;
696      } else if(rowActivity[i]>rowUpper[i]) {
697        if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
698          feasible = false;
699      }
700    }
701    for (int i=0; i<numberIntegers; i++) {
702      int iColumn = integerVariable[i];
703      double value=newSolution[iColumn];
704      if (fabs(floor(value+0.5)-value)>integerTolerance) {
705        feasible = false;
706        break;
707      }
708    }
709    if (feasible) {
710      // new solution
711      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
712      solutionValue = newSolutionValue;
713      //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
714      returnCode=1;
715    } else {
716      // Can easily happen
717      //printf("Debug CbcHeuristicDive giving bad solution\n");
718    }
719  }
720
721#ifdef DIVE_DEBUG
722  std::cout<<"nRoundInfeasible = "<<nRoundInfeasible
723           <<", nRoundFeasible = "<<nRoundFeasible
724           <<", returnCode = "<<returnCode
725           <<", reasonToStop = "<<reasonToStop
726           <<", simplexIts = "<<numberSimplexIterations
727           <<", iterations = "<<iteration<<std::endl;
728#endif
729
730  delete [] newSolution;
731  delete [] columnFixed;
732  delete [] originalBound;
733  delete [] fixedAtLowerBound;
734  delete [] candidate;
735  delete [] rowActivity;
736  delete [] random;
737  delete [] downArray_;
738  downArray_=NULL;
739  delete [] upArray_;
740  upArray_=NULL;
741  delete solver;
742  return returnCode;
743}
744
745// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
746void 
747CbcHeuristicDive::validate() 
748{
749  if (model_&&(when()%100)<10) {
750    if (model_->numberIntegers()!=
751        model_->numberObjects()&&(model_->numberObjects()||
752                                  (model_->specialOptions()&1024)==0)) {
753      int numberOdd=0;
754      for (int i=0;i<model_->numberObjects();i++) {
755        if (!model_->object(i)->canDoHeuristics()) 
756          numberOdd++;
757      }
758      if (numberOdd)
759        setWhen(0);
760    }
761  }
762
763  int numberIntegers = model_->numberIntegers();
764  const int * integerVariable = model_->integerVariable();
765  delete [] downLocks_;
766  delete [] upLocks_;
767  downLocks_ = new unsigned short [numberIntegers];
768  upLocks_ = new unsigned short [numberIntegers];
769  // Column copy
770  const double * element = matrix_.getElements();
771  const int * row = matrix_.getIndices();
772  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
773  const int * columnLength = matrix_.getVectorLengths();
774  const double * rowLower = model_->solver()->getRowLower();
775  const double * rowUpper = model_->solver()->getRowUpper();
776  for (int i=0;i<numberIntegers;i++) {
777    int iColumn = integerVariable[i];
778    int down=0;
779    int up=0;
780    if (columnLength[iColumn]>65535) {
781      setWhen(0);
782      break; // unlikely to work
783    }
784    for (CoinBigIndex j=columnStart[iColumn];
785         j<columnStart[iColumn]+columnLength[iColumn];j++) {
786      int iRow=row[j];
787      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
788        up++;
789        down++;
790      } else if (element[j]>0.0) {
791        if (rowUpper[iRow]<1.0e20)
792          up++;
793        else
794          down++;
795      } else {
796        if (rowLower[iRow]>-1.0e20)
797          up++;
798        else
799          down++;
800      }
801    }
802    downLocks_[i] = static_cast<unsigned short> (down);
803    upLocks_[i] = static_cast<unsigned short> (up);
804  }
805
806#ifdef DIVE_FIX_BINARY_VARIABLES
807  selectBinaryVariables();
808#endif
809}
810
811// Select candidate binary variables for fixing
812void
813CbcHeuristicDive::selectBinaryVariables()
814{
815  // Row copy
816  const double * elementByRow = matrixByRow_.getElements();
817  const int * column = matrixByRow_.getIndices();
818  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
819  const int * rowLength = matrixByRow_.getVectorLengths();
820
821  const int numberRows = matrixByRow_.getNumRows();
822  const int numberCols = matrixByRow_.getNumCols();
823
824  const double * lower = model_->solver()->getColLower();
825  const double * upper = model_->solver()->getColUpper();
826  const double * rowLower = model_->solver()->getRowLower();
827  const double * rowUpper = model_->solver()->getRowUpper();
828
829  //  const char * integerType = model_->integerType();
830 
831
832  //  const int numberIntegers = model_->numberIntegers();
833  //  const int * integerVariable = model_->integerVariable();
834  const double * objective = model_->solver()->getObjCoefficients();
835
836  // vector to store the row number of variable bound rows
837  int* rowIndexes = new int [numberCols];
838  memset(rowIndexes, -1, numberCols*sizeof(int));
839
840  for(int i=0; i<numberRows; i++) {
841    int positiveBinary = -1;
842    int negativeBinary = -1;
843    int nPositiveOther = 0;
844    int nNegativeOther = 0;
845    for (int k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
846      int iColumn = column[k];
847      if(model_->solver()->isInteger(iColumn) &&
848         lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
849         objective[iColumn] == 0.0 &&
850         elementByRow[k] > 0.0 &&
851         positiveBinary < 0)
852        positiveBinary = iColumn;
853      else if(model_->solver()->isInteger(iColumn) &&
854              lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
855              objective[iColumn] == 0.0 &&
856              elementByRow[k] < 0.0 &&
857              negativeBinary < 0)
858        negativeBinary = iColumn;
859      else if((elementByRow[k]> 0.0 &&
860               lower[iColumn] >= 0.0) ||
861              (elementByRow[k] < 0.0 &&
862               upper[iColumn] <= 0.0))
863        nPositiveOther++;
864      else if((elementByRow[k]> 0.0 &&
865               lower[iColumn] <= 0.0) ||
866              (elementByRow[k] < 0.0 &&
867               upper[iColumn] >= 0.0))
868        nNegativeOther++;
869      if(nPositiveOther > 0 && nNegativeOther > 0)
870        break;
871    }
872    int binVar = -1;
873    if(positiveBinary >= 0 &&
874       (negativeBinary >= 0 || nNegativeOther > 0) &&
875       nPositiveOther == 0 &&
876       rowLower[i] == 0.0 &&
877       rowUpper[i] > 0.0)
878      binVar = positiveBinary;
879    else if(negativeBinary >= 0 &&
880            (positiveBinary >= 0 || nPositiveOther > 0) &&
881            nNegativeOther == 0 &&
882            rowLower[i] < 0.0 &&
883            rowUpper[i] == 0.0)
884      binVar = negativeBinary;
885    if(binVar >= 0) {
886      if(rowIndexes[binVar] == -1)
887        rowIndexes[binVar] = i;
888      else if(rowIndexes[binVar] >= 0)
889        rowIndexes[binVar] = -2;
890    }
891  }
892
893  for(int j=0; j<numberCols; j++) {
894    if(rowIndexes[j] >= 0) {
895      binVarIndex_.push_back(j);
896      vbRowIndex_.push_back(rowIndexes[j]);
897    }
898  }
899
900#ifdef DIVE_DEBUG
901  std::cout<<"number vub Binary = "<<binVarIndex_.size()<<std::endl;
902#endif
903
904  delete [] rowIndexes;
905   
906}
907
908/*
909  Perform reduced cost fixing on integer variables.
910
911  The variables in question are already nonbasic at bound. We're just nailing
912  down the current situation.
913*/
914
915int CbcHeuristicDive::reducedCostFix (OsiSolverInterface* solver)
916
917{
918  //return 0; // temp
919#if 1
920  if(!model_->solverCharacteristics()->reducedCostsAccurate())
921    return 0; //NLP
922#endif
923  double cutoff = model_->getCutoff() ;
924  if (cutoff>1.0e20)
925    return 0;
926#ifdef DIVE_DEBUG
927  std::cout<<"cutoff = "<<cutoff<<std::endl;
928#endif
929  double direction = solver->getObjSense() ;
930  double gap = cutoff - solver->getObjValue()*direction ;
931  gap *= 0.5; // Fix more
932  double tolerance;
933  solver->getDblParam(OsiDualTolerance,tolerance) ;
934  if (gap<=0.0)
935    gap = tolerance; //return 0;
936  gap += 100.0*tolerance;
937  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
938
939  const double *lower = solver->getColLower() ;
940  const double *upper = solver->getColUpper() ;
941  const double *solution = solver->getColSolution() ;
942  const double *reducedCost = solver->getReducedCost() ;
943
944  int numberIntegers = model_->numberIntegers();
945  const int * integerVariable = model_->integerVariable();
946
947  int numberFixed = 0 ;
948
949# ifdef COIN_HAS_CLP
950  OsiClpSolverInterface * clpSolver
951    = dynamic_cast<OsiClpSolverInterface *> (solver);
952  ClpSimplex * clpSimplex=NULL;
953  if (clpSolver) 
954    clpSimplex = clpSolver->getModelPtr();
955# endif
956  for (int i = 0 ; i < numberIntegers ; i++) {
957    int iColumn = integerVariable[i] ;
958    double djValue = direction*reducedCost[iColumn] ;
959    if (upper[iColumn]-lower[iColumn] > integerTolerance) {
960      if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap) {
961#ifdef COIN_HAS_CLP
962        // may just have been fixed before
963        if (clpSimplex) {
964          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
965#ifdef COIN_DEVELOP
966            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
967                   iColumn,clpSimplex->getColumnStatus(iColumn),
968                   djValue,gap,lower[iColumn],upper[iColumn]);
969#endif
970          } else {         
971            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
972                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
973          }
974        }
975#endif
976        solver->setColUpper(iColumn,lower[iColumn]) ;
977        numberFixed++ ;
978      } else if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap) {
979#ifdef COIN_HAS_CLP
980        // may just have been fixed before
981        if (clpSimplex) {
982          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
983#ifdef COIN_DEVELOP
984            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
985                   iColumn,clpSimplex->getColumnStatus(iColumn),
986                   djValue,gap,lower[iColumn],upper[iColumn]);
987#endif
988          } else {         
989            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
990                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
991          }
992        }
993#endif
994        solver->setColLower(iColumn,upper[iColumn]) ;
995        numberFixed++ ;
996      }
997    }
998  }
999  return numberFixed; 
1000}
1001// Fix other variables at bounds
1002int 
1003CbcHeuristicDive::fixOtherVariables(OsiSolverInterface * solver,
1004                                    const double * solution,
1005                                    PseudoReducedCost * candidate,
1006                                    const double * random)
1007{
1008  const double * lower = solver->getColLower();
1009  const double * upper = solver->getColUpper();
1010  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1011  double primalTolerance;
1012  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
1013
1014  int numberIntegers = model_->numberIntegers();
1015  const int * integerVariable = model_->integerVariable();
1016  const double* reducedCost = solver->getReducedCost();
1017  // fix other integer variables that are at their bounds
1018  int cnt=0;
1019#ifdef GAP
1020  double direction = solver->getObjSense(); // 1 for min, -1 for max
1021  double gap=1.0e30;
1022#endif
1023#ifdef GAP
1024  double cutoff = model_->getCutoff() ; 
1025  if (cutoff<1.0e20&&false) {
1026    double direction = solver->getObjSense() ;
1027    gap = cutoff - solver->getObjValue()*direction ;
1028    gap *= 0.1; // Fix more if plausible
1029    double tolerance;
1030    solver->getDblParam(OsiDualTolerance,tolerance) ;
1031    if (gap<=0.0)
1032      gap = tolerance;
1033    gap += 100.0*tolerance;
1034  }
1035  int nOverGap=0;
1036#endif
1037  int numberFree=0;
1038  int numberFixedAlready=0;
1039  for (int i=0; i<numberIntegers; i++) {
1040    int iColumn = integerVariable[i];
1041    if (upper[iColumn]>lower[iColumn]) {
1042      numberFree++;
1043      double value=solution[iColumn];
1044      if(fabs(floor(value+0.5)-value)<=integerTolerance) {
1045        candidate[cnt].var = iColumn;
1046        candidate[cnt++].pseudoRedCost = 
1047          fabs(reducedCost[iColumn]*random[i]);
1048#ifdef GAP
1049        if (fabs(reducedCost[iColumn])>gap)
1050          nOverGap++;
1051#endif
1052      }
1053    } else {
1054      numberFixedAlready++;
1055    }
1056  }
1057#ifdef GAP
1058  int nLeft = maxNumberToFix-numberFixedAlready;
1059#ifdef CLP_INVESTIGATE
1060  printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
1061         cutoff,solver->getObjValue(),nOverGap,numberFree,
1062         numberFixedAlready);
1063#endif
1064  if (nOverGap>nLeft&&true) {
1065    nOverGap = CoinMin(nOverGap,nLeft+maxNumberToFix/2);
1066    maxNumberToFix += nOverGap-nLeft;
1067  }
1068#else
1069#ifdef CLP_INVESTIGATE
1070  printf("cutoff %g obj %g - %d free, %d fixed\n",
1071         model_->getCutoff(),solver->getObjValue(),numberFree,
1072         numberFixedAlready);
1073#endif
1074#endif
1075  return cnt;
1076}
Note: See TracBrowser for help on using the repository browser.