Ignore:
Timestamp:
Apr 10, 2008 1:55:10 PM (13 years ago)
Author:
ladanyi
Message:

Incorporated changes from branches/heur

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcHeuristicDiveGuided.cpp

    r871 r912  
    88#include "CbcHeuristicDiveGuided.hpp"
    99#include "CbcStrategy.hpp"
    10 #include  "CoinTime.hpp"
    1110
    1211// Default Constructor
    1312CbcHeuristicDiveGuided::CbcHeuristicDiveGuided()
    14   :CbcHeuristic()
     13  :CbcHeuristicDive()
    1514{
    16   // matrix and row copy will automatically be empty
    17   downLocks_ =NULL;
    18   upLocks_ = NULL;
    19   percentageToFix_ = 0.2;
    20   maxIterations_ = 100;
    21   maxTime_ = 60;
    2215}
    2316
    2417// Constructor from model
    2518CbcHeuristicDiveGuided::CbcHeuristicDiveGuided(CbcModel & model)
    26   :CbcHeuristic(model)
     19  :CbcHeuristicDive(model)
    2720{
    28   downLocks_ =NULL;
    29   upLocks_ = NULL;
    30   // Get a copy of original matrix
    31   assert(model.solver());
    32   // model may have empty matrix - wait until setModel
    33   const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
    34   if (matrix) {
    35     matrix_ = *matrix;
    36   }
    37   percentageToFix_ = 0.2;
    38   maxIterations_ = 100;
    39   maxTime_ = 60;
    4021}
    4122
     
    4324CbcHeuristicDiveGuided::~CbcHeuristicDiveGuided ()
    4425{
    45   delete [] downLocks_;
    46   delete [] upLocks_;
    4726}
    4827
     
    6241  fprintf(fp,"3  CbcHeuristicDiveGuided heuristicDiveGuided(*cbcModel);\n");
    6342  CbcHeuristic::generateCpp(fp,"heuristicDiveGuided");
    64   if (percentageToFix_!=other.percentageToFix_)
    65     fprintf(fp,"3  heuristicDiveGuided.setPercentageToFix(%.2f);\n",percentageToFix_);
    66   else
    67     fprintf(fp,"4  heuristicDiveGuided.setPercentageToFix(%.2f);\n",percentageToFix_);
    68   if (maxIterations_!=other.maxIterations_)
    69     fprintf(fp,"3  heuristicDiveGuided.setMaxIterations(%d);\n",maxIterations_);
    70   else
    71     fprintf(fp,"4  heuristicDiveGuided.setMaxIterations(%d);\n",maxIterations_);
    72   if (maxTime_!=other.maxTime_)
    73     fprintf(fp,"3  heuristicDiveGuided.setMaxTime(%.2f);\n",maxTime_);
    74   else
    75     fprintf(fp,"4  heuristicDiveGuided.setMaxTime(%.2f);\n",maxTime_);
    7643  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveGuided);\n");
    7744}
     
    8047CbcHeuristicDiveGuided::CbcHeuristicDiveGuided(const CbcHeuristicDiveGuided & rhs)
    8148:
    82   CbcHeuristic(rhs),
    83   matrix_(rhs.matrix_),
    84   percentageToFix_(rhs.percentageToFix_),
    85   maxIterations_(rhs.maxIterations_),
    86   maxTime_(rhs.maxTime_)
     49  CbcHeuristicDive(rhs)
    8750{
    88   int numberIntegers = model_->numberIntegers();
    89   if (rhs.downLocks_) {
    90     int numberIntegers = model_->numberIntegers();
    91     downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    92     upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    93   } else {
    94     downLocks_ = NULL;
    95     upLocks_ = NULL;
    96   }
    9751}
    9852
     
    10256{
    10357  if (this!=&rhs) {
    104     CbcHeuristic::operator=(rhs);
    105     matrix_ = rhs.matrix_;
    106     percentageToFix_ = rhs.percentageToFix_;
    107     maxIterations_ = rhs.maxIterations_;
    108     maxTime_ = rhs.maxTime_;
    109     delete [] downLocks_;
    110     delete [] upLocks_;
    111     if (rhs.downLocks_) {
    112       int numberIntegers = model_->numberIntegers();
    113       downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    114       upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    115     } else {
    116       downLocks_ = NULL;
    117       upLocks_ = NULL;
    118     }
     58    CbcHeuristicDive::operator=(rhs);
    11959  }
    12060  return *this;
    12161}
    12262
    123 // Resets stuff if model changes
    124 void
    125 CbcHeuristicDiveGuided::resetModel(CbcModel * model)
     63bool
     64CbcHeuristicDiveGuided::canHeuristicRun()
    12665{
    127   model_=model;
    128   // Get a copy of original matrix
    129   assert(model_->solver());
    130   // model may have empty matrix - wait until setModel
    131   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    132   if (matrix) {
    133     matrix_ = *matrix;
    134     validate();
    135   }
     66  double* bestIntegerSolution = model_->bestSolution();
     67  if(bestIntegerSolution == NULL)
     68    return false; // no integer solution available. Switch off heuristic
     69
     70  return CbcHeuristicDive::canHeuristicRun();
    13671}
    13772
    138 // See if dive guided will give better solution
    139 // Sets value of solution
    140 // Returns 1 if solution, 0 if not
    141 int
    142 CbcHeuristicDiveGuided::solution(double & solutionValue,
    143                                  double * betterSolution)
     73void
     74CbcHeuristicDiveGuided::selectVariableToBranch(OsiSolverInterface* solver,
     75                                                   const double* newSolution,
     76                                                   int& bestColumn,
     77                                                   int& bestRound)
    14478{
     79  double* bestIntegerSolution = model_->bestSolution();
    14580
    146   // See if to do
    147   if (!when()||(when()%10==1&&model_->phase()!=1)||
    148       (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
    149     return 0; // switched off
    150 
    151   double* bestIntegerSolution = model_->bestSolution();
    152   if(bestIntegerSolution == NULL)
    153     return 0; // no integer solution available. Switch off heuristic
    154 
    155   double time1 = CoinCpuTime();
    156 
    157   OsiSolverInterface * solver = model_->solver()->clone();
    158   const double * lower = solver->getColLower();
    159   const double * upper = solver->getColUpper();
    160   const double * rowLower = solver->getRowLower();
    161   const double * rowUpper = solver->getRowUpper();
    162   const double * solution = solver->getColSolution();
    163   const double * objective = solver->getObjCoefficients();
    164   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    165   double primalTolerance;
    166   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    167 
    168   int numberRows = matrix_.getNumRows();
    169   assert (numberRows<=solver->getNumRows());
    17081  int numberIntegers = model_->numberIntegers();
    17182  const int * integerVariable = model_->integerVariable();
    172   double direction = solver->getObjSense(); // 1 for min, -1 for max
    173   double newSolutionValue = direction*solver->getObjValue();
    174   int returnCode = 0;
    175   // Column copy
    176   const double * element = matrix_.getElements();
    177   const int * row = matrix_.getIndices();
    178   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    179   const int * columnLength = matrix_.getVectorLengths();
     83  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    18084
    181   // Get solution array for heuristic solution
    182   int numberColumns = solver->getNumCols();
    183   double * newSolution = new double [numberColumns];
    184   memcpy(newSolution,solution,numberColumns*sizeof(double));
    185 
    186   // vectors to store the latest variables fixed at their bounds
    187   int* columnFixed = new int [numberIntegers];
    188   double* originalBound = new double [numberIntegers];
    189   bool * fixedAtLowerBound = new bool [numberIntegers];
    190 
    191   const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
    192 
    193   // count how many fractional variables
    194   int numberFractionalVariables = 0;
     85  bestColumn = -1;
     86  bestRound = -1; // -1 rounds down, +1 rounds up
     87  double bestFraction = DBL_MAX;
    19588  for (int i=0; i<numberIntegers; i++) {
    19689    int iColumn = integerVariable[i];
    19790    double value=newSolution[iColumn];
     91    double fraction=value-floor(value);
     92    int round = 0;
    19893    if (fabs(floor(value+0.5)-value)>integerTolerance) {
    199       numberFractionalVariables++;
    200     }
    201   }
     94      if(downLocks_[i]>0&&upLocks_[i]>0) {
     95        // the variable cannot be rounded
     96        if(value >= bestIntegerSolution[iColumn])
     97          round = -1;
     98        else {
     99          round = 1;
     100          fraction = 1.0 - fraction;
     101        }
     102       
     103        // if variable is not binary, penalize it
     104        if(!solver->isBinary(iColumn))
     105          fraction *= 1000.0;
    202106
    203   int iteration = 0;
    204   while(numberFractionalVariables) {
    205     iteration++;
    206 
    207     // select a fractional variable to bound
    208     double bestFraction = DBL_MAX;
    209     int bestColumn = -1;
    210     int bestRound = -1; // -1 rounds down, +1 rounds up
    211     int numberAtBoundFixed = 0;
    212     bool canRoundSolution = true;
    213     for (int i=0; i<numberIntegers; i++) {
    214       int iColumn = integerVariable[i];
    215       double value=newSolution[iColumn];
    216       double fraction=value-floor(value);
    217       int round = 0;
    218       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    219         if(downLocks_[i]>0&&upLocks_[i]>0) {
    220           // the variable cannot be rounded
    221           canRoundSolution = false;
    222           if(value >= bestIntegerSolution[iColumn])
    223             round = -1;
    224           else {
    225             round = 1;
    226             fraction = 1.0 - fraction;
    227           }
    228 
    229           // if variable is not binary, penalize it
    230           if(!solver->isBinary(iColumn))
    231             fraction *= 1000.0;
    232 
    233           if(fraction < bestFraction) {
    234             bestColumn = iColumn;
    235             bestFraction = fraction;
    236             bestRound = round;
    237           }
    238         }
    239       }
    240       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    241         // fix the variable at one of its bounds
    242         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    243             lower[iColumn] != upper[iColumn]) {
    244           columnFixed[numberAtBoundFixed] = iColumn;
    245           originalBound[numberAtBoundFixed] = upper[iColumn];
    246           fixedAtLowerBound[numberAtBoundFixed] = true;
    247           solver->setColUpper(iColumn, lower[iColumn]);
    248           numberAtBoundFixed++;
    249         }
    250         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    251             lower[iColumn] != upper[iColumn]) {
    252           columnFixed[numberAtBoundFixed] = iColumn;
    253           originalBound[numberAtBoundFixed] = lower[iColumn];
    254           fixedAtLowerBound[numberAtBoundFixed] = false;
    255           solver->setColLower(iColumn, upper[iColumn]);
    256           numberAtBoundFixed++;
     107        if(fraction < bestFraction) {
     108          bestColumn = iColumn;
     109          bestFraction = fraction;
     110          bestRound = round;
    257111        }
    258112      }
    259113    }
    260 
    261     if(canRoundSolution) {
    262       // Round all the fractional variables
    263       for (int i=0; i<numberIntegers; i++) {
    264         int iColumn = integerVariable[i];
    265         double value=newSolution[iColumn];
    266         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    267           if(downLocks_[i]==0 || upLocks_[i]==0) {
    268             if(downLocks_[i]==0 && upLocks_[i]==0) {
    269               if(direction * objective[iColumn] >= 0.0)
    270                 newSolution[iColumn] = floor(value);
    271               else
    272                 newSolution[iColumn] = ceil(value);
    273             }
    274             else if(downLocks_[i]==0)
    275               newSolution[iColumn] = floor(value);
    276             else
    277               newSolution[iColumn] = ceil(value);
    278           }
    279           else
    280             break;
    281         }
    282       }
    283       break;
    284     }
    285 
    286 
    287     double originalBoundBestColumn;
    288     if(bestColumn >= 0) {
    289       if(bestRound < 0) {
    290         originalBoundBestColumn = upper[bestColumn];
    291         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    292       }
    293       else {
    294         originalBoundBestColumn = lower[bestColumn];
    295         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    296       }
    297     } else {
    298       break;
    299     }
    300     int originalBestRound = bestRound;
    301     while (1) {
    302 
    303       solver->resolve();
    304 
    305       if(!solver->isProvenOptimal()) {
    306         if(numberAtBoundFixed > 0) {
    307           // Remove the bound fix for variables that were at bounds
    308           for(int i=0; i<numberAtBoundFixed; i++) {
    309             int iColFixed = columnFixed[i];
    310             if(fixedAtLowerBound[i])
    311               solver->setColUpper(iColFixed, originalBound[i]);
    312             else
    313               solver->setColLower(iColFixed, originalBound[i]);
    314           }
    315           numberAtBoundFixed = 0;
    316         }
    317         else if(bestRound == originalBestRound) {
    318           bestRound *= (-1);
    319           if(bestRound < 0) {
    320             solver->setColLower(bestColumn, originalBoundBestColumn);
    321             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    322           }
    323           else {
    324             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    325             solver->setColUpper(bestColumn, originalBoundBestColumn);
    326           }
    327         }
    328         else
    329           break;
    330       }
    331       else
    332         break;
    333     }
    334 
    335     if(!solver->isProvenOptimal())
    336       break;
    337 
    338     if(iteration > maxIterations_) {
    339       break;
    340     }
    341 
    342     if(CoinCpuTime()-time1 > maxTime_) {
    343       break;
    344     }
    345 
    346     memcpy(newSolution,solution,numberColumns*sizeof(double));
    347     numberFractionalVariables = 0;
    348     for (int i=0; i<numberIntegers; i++) {
    349       int iColumn = integerVariable[i];
    350       double value=newSolution[iColumn];
    351       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    352         numberFractionalVariables++;
    353       }
    354     }
    355 
    356   }
    357 
    358 
    359   double * rowActivity = new double[numberRows];
    360   memset(rowActivity,0,numberRows*sizeof(double));
    361 
    362   // re-compute new solution value
    363   double objOffset=0.0;
    364   solver->getDblParam(OsiObjOffset,objOffset);
    365   newSolutionValue = -objOffset;
    366   for (int i=0 ; i<numberColumns ; i++ )
    367     newSolutionValue += objective[i]*newSolution[i];
    368   newSolutionValue *= direction;
    369     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    370   if (newSolutionValue<solutionValue) {
    371     // paranoid check
    372     memset(rowActivity,0,numberRows*sizeof(double));
    373     for (int i=0;i<numberColumns;i++) {
    374       int j;
    375       double value = newSolution[i];
    376       if (value) {
    377         for (j=columnStart[i];
    378              j<columnStart[i]+columnLength[i];j++) {
    379           int iRow=row[j];
    380           rowActivity[iRow] += value*element[j];
    381         }
    382       }
    383     }
    384     // check was approximately feasible
    385     bool feasible=true;
    386     for (int i=0;i<numberRows;i++) {
    387       if(rowActivity[i]<rowLower[i]) {
    388         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    389           feasible = false;
    390       } else if(rowActivity[i]>rowUpper[i]) {
    391         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    392           feasible = false;
    393       }
    394     }
    395     for (int i=0; i<numberIntegers; i++) {
    396       int iColumn = integerVariable[i];
    397       double value=newSolution[iColumn];
    398       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    399         feasible = false;
    400         break;
    401       }
    402     }
    403     if (feasible) {
    404       // new solution
    405       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    406       solutionValue = newSolutionValue;
    407       //printf("** Solution of %g found by CbcHeuristicDiveGuided\n",newSolutionValue);
    408       returnCode=1;
    409     } else {
    410       // Can easily happen
    411       //printf("Debug CbcHeuristicDiveGuided giving bad solution\n");
    412     }
    413   }
    414 
    415   delete [] newSolution;
    416   delete [] columnFixed;
    417   delete [] originalBound;
    418   delete [] fixedAtLowerBound;
    419   delete [] rowActivity;
    420   delete solver;
    421   return returnCode;
    422 }
    423 
    424 // update model
    425 void CbcHeuristicDiveGuided::setModel(CbcModel * model)
    426 {
    427   model_ = model;
    428   // Get a copy of original matrix (and by row for rounding);
    429   assert(model_->solver());
    430   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    431   if (matrix) {
    432     matrix_ = *matrix;
    433     // make sure model okay for heuristic
    434     validate();
    435114  }
    436115}
    437 
    438 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    439 void
    440 CbcHeuristicDiveGuided::validate()
    441 {
    442   if (model_&&when()<10) {
    443     if (model_->numberIntegers()!=
    444         model_->numberObjects())
    445       setWhen(0);
    446   }
    447 
    448   int numberIntegers = model_->numberIntegers();
    449   const int * integerVariable = model_->integerVariable();
    450   delete [] downLocks_;
    451   delete [] upLocks_;
    452   downLocks_ = new unsigned short [numberIntegers];
    453   upLocks_ = new unsigned short [numberIntegers];
    454   // Column copy
    455   const double * element = matrix_.getElements();
    456   const int * row = matrix_.getIndices();
    457   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    458   const int * columnLength = matrix_.getVectorLengths();
    459   const double * rowLower = model_->solver()->getRowLower();
    460   const double * rowUpper = model_->solver()->getRowUpper();
    461   for (int i=0;i<numberIntegers;i++) {
    462     int iColumn = integerVariable[i];
    463     int down=0;
    464     int up=0;
    465     if (columnLength[iColumn]>65535) {
    466       setWhen(0);
    467       break; // unlikely to work
    468     }
    469     for (CoinBigIndex j=columnStart[iColumn];
    470          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    471       int iRow=row[j];
    472       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    473         up++;
    474         down++;
    475       } else if (element[j]>0.0) {
    476         if (rowUpper[iRow]<1.0e20)
    477           up++;
    478         else
    479           down++;
    480       } else {
    481         if (rowLower[iRow]>-1.0e20)
    482           up++;
    483         else
    484           down++;
    485       }
    486     }
    487     downLocks_[i] = (unsigned short) down;
    488     upLocks_[i] = (unsigned short) up;
    489   }
    490 }
Note: See TracChangeset for help on using the changeset viewer.