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/CbcHeuristicDiveFractional.cpp

    r871 r912  
    88#include "CbcHeuristicDiveFractional.hpp"
    99#include "CbcStrategy.hpp"
    10 #include  "CoinTime.hpp"
    1110
    1211// Default Constructor
    1312CbcHeuristicDiveFractional::CbcHeuristicDiveFractional()
    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
    2518CbcHeuristicDiveFractional::CbcHeuristicDiveFractional(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
     
    4324CbcHeuristicDiveFractional::~CbcHeuristicDiveFractional ()
    4425{
    45   delete [] downLocks_;
    46   delete [] upLocks_;
    4726}
    4827
     
    6241  fprintf(fp,"3  CbcHeuristicDiveFractional heuristicDiveFractional(*cbcModel);\n");
    6342  CbcHeuristic::generateCpp(fp,"heuristicDiveFractional");
    64   if (percentageToFix_!=other.percentageToFix_)
    65     fprintf(fp,"3  heuristicDiveFractional.setPercentageToFix(%.f);\n",percentageToFix_);
    66   else
    67     fprintf(fp,"4  heuristicDiveFractional.setPercentageToFix(%.f);\n",percentageToFix_);
    68   if (maxIterations_!=other.maxIterations_)
    69     fprintf(fp,"3  heuristicDiveFractional.setMaxIterations(%d);\n",maxIterations_);
    70   else
    71     fprintf(fp,"4  heuristicDiveFractional.setMaxIterations(%d);\n",maxIterations_);
    72   if (maxTime_!=other.maxTime_)
    73     fprintf(fp,"3  heuristicDiveFractional.setMaxTime(%.2f);\n",maxTime_);
    74   else
    75     fprintf(fp,"4  heuristicDiveFractional.setMaxTime(%.2f);\n",maxTime_);
    7643  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveFractional);\n");
    7744}
     
    8047CbcHeuristicDiveFractional::CbcHeuristicDiveFractional(const CbcHeuristicDiveFractional & 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 CbcHeuristicDiveFractional::resetModel(CbcModel * model)
     63void
     64CbcHeuristicDiveFractional::selectVariableToBranch(OsiSolverInterface* solver,
     65                                                   const double* newSolution,
     66                                                   int& bestColumn,
     67                                                   int& bestRound)
    12668{
    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   }
    136 }
    137 
    138 // See if dive fractional will give better solution
    139 // Sets value of solution
    140 // Returns 1 if solution, 0 if not
    141 int
    142 CbcHeuristicDiveFractional::solution(double & solutionValue,
    143                                      double * betterSolution)
    144 {
    145 
    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 time1 = CoinCpuTime();
    152 
    153   OsiSolverInterface * solver = model_->solver()->clone();
    154   const double * lower = solver->getColLower();
    155   const double * upper = solver->getColUpper();
    156   const double * rowLower = solver->getRowLower();
    157   const double * rowUpper = solver->getRowUpper();
    158   const double * solution = solver->getColSolution();
    159   const double * objective = solver->getObjCoefficients();
    160   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    161   double primalTolerance;
    162   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    163 
    164   int numberRows = matrix_.getNumRows();
    165   assert (numberRows<=solver->getNumRows());
    16669  int numberIntegers = model_->numberIntegers();
    16770  const int * integerVariable = model_->integerVariable();
    168   double direction = solver->getObjSense(); // 1 for min, -1 for max
    169   double newSolutionValue = direction*solver->getObjValue();
    170   int returnCode = 0;
    171   // Column copy
    172   const double * element = matrix_.getElements();
    173   const int * row = matrix_.getIndices();
    174   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    175   const int * columnLength = matrix_.getVectorLengths();
     71  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    17672
    177   // Get solution array for heuristic solution
    178   int numberColumns = solver->getNumCols();
    179   double * newSolution = new double [numberColumns];
    180   memcpy(newSolution,solution,numberColumns*sizeof(double));
    181 
    182   // vectors to store the latest variables fixed at their bounds
    183   int* columnFixed = new int [numberIntegers];
    184   double* originalBound = new double [numberIntegers];
    185   bool * fixedAtLowerBound = new bool [numberIntegers];
    186 
    187   const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
    188 
    189   // count how many fractional variables
    190   int numberFractionalVariables = 0;
     73  bestColumn = -1;
     74  bestRound = -1; // -1 rounds down, +1 rounds up
     75  double bestFraction = DBL_MAX;
    19176  for (int i=0; i<numberIntegers; i++) {
    19277    int iColumn = integerVariable[i];
    19378    double value=newSolution[iColumn];
     79    double fraction=value-floor(value);
     80    int round = 0;
    19481    if (fabs(floor(value+0.5)-value)>integerTolerance) {
    195       numberFractionalVariables++;
    196     }
    197   }
    198 
    199   int iteration = 0;
    200   while(numberFractionalVariables) {
    201     iteration++;
    202 
    203     // select a fractional variable to bound
    204     double bestFraction = DBL_MAX;
    205     int bestColumn = -1;
    206     int bestRound = -1; // -1 rounds down, +1 rounds up
    207     int numberAtBoundFixed = 0;
    208     bool canRoundSolution = true;
    209     for (int i=0; i<numberIntegers; i++) {
    210       int iColumn = integerVariable[i];
    211       double value=newSolution[iColumn];
    212       double fraction=value-floor(value);
    213       int round = 0;
    214       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    215         if(downLocks_[i]>0&&upLocks_[i]>0) {
     82      if(downLocks_[i]>0&&upLocks_[i]>0) {
    21683          // the variable cannot be rounded
    217           canRoundSolution = false;
    218           if(fraction < 0.5)
    219             round = -1;
    220           else {
    221             round = 1;
    222             fraction = 1.0 - fraction;
    223           }
    224 
    225           // if variable is not binary, penalize it
    226           if(!solver->isBinary(iColumn))
     84        if(fraction < 0.5)
     85          round = -1;
     86        else {
     87          round = 1;
     88          fraction = 1.0 - fraction;
     89        }
     90       
     91        // if variable is not binary, penalize it
     92        if(!solver->isBinary(iColumn))
    22793            fraction *= 1000.0;
    228 
    229           if(fraction < bestFraction) {
    230             bestColumn = iColumn;
    231             bestFraction = fraction;
    232             bestRound = round;
    233           }
    234         }
    235       }
    236       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    237         // fix the variable at one of its bounds
    238         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    239             lower[iColumn] != upper[iColumn]) {
    240           columnFixed[numberAtBoundFixed] = iColumn;
    241           originalBound[numberAtBoundFixed] = upper[iColumn];
    242           fixedAtLowerBound[numberAtBoundFixed] = true;
    243           solver->setColUpper(iColumn, lower[iColumn]);
    244           numberAtBoundFixed++;
    245         }
    246         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    247             lower[iColumn] != upper[iColumn]) {
    248           columnFixed[numberAtBoundFixed] = iColumn;
    249           originalBound[numberAtBoundFixed] = lower[iColumn];
    250           fixedAtLowerBound[numberAtBoundFixed] = false;
    251           solver->setColLower(iColumn, upper[iColumn]);
    252           numberAtBoundFixed++;
     94       
     95        if(fraction < bestFraction) {
     96          bestColumn = iColumn;
     97          bestFraction = fraction;
     98          bestRound = round;
    25399        }
    254100      }
    255101    }
    256 
    257     if(canRoundSolution) {
    258       // Round all the fractional variables
    259       for (int i=0; i<numberIntegers; i++) {
    260         int iColumn = integerVariable[i];
    261         double value=newSolution[iColumn];
    262         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    263           if(downLocks_[i]==0 || upLocks_[i]==0) {
    264             if(downLocks_[i]==0 && upLocks_[i]==0) {
    265               if(direction * objective[iColumn] >= 0.0)
    266                 newSolution[iColumn] = floor(value);
    267               else
    268                 newSolution[iColumn] = ceil(value);
    269             }
    270             else if(downLocks_[i]==0)
    271               newSolution[iColumn] = floor(value);
    272             else
    273               newSolution[iColumn] = ceil(value);
    274           }
    275           else
    276             break;
    277         }
    278       }
    279       break;
    280     }
    281 
    282 
    283     double originalBoundBestColumn;
    284     if(bestColumn >= 0) {
    285       if(bestRound < 0) {
    286         originalBoundBestColumn = upper[bestColumn];
    287         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    288       }
    289       else {
    290         originalBoundBestColumn = lower[bestColumn];
    291         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    292       }
    293     } else {
    294       break;
    295     }
    296     int originalBestRound = bestRound;
    297     while (1) {
    298 
    299       solver->resolve();
    300 
    301       if(!solver->isProvenOptimal()) {
    302         if(numberAtBoundFixed > 0) {
    303           // Remove the bound fix for variables that were at bounds
    304           for(int i=0; i<numberAtBoundFixed; i++) {
    305             int iColFixed = columnFixed[i];
    306             if(fixedAtLowerBound[i])
    307               solver->setColUpper(iColFixed, originalBound[i]);
    308             else
    309               solver->setColLower(iColFixed, originalBound[i]);
    310           }
    311           numberAtBoundFixed = 0;
    312         }
    313         else if(bestRound == originalBestRound) {
    314           bestRound *= (-1);
    315           if(bestRound < 0) {
    316             solver->setColLower(bestColumn, originalBoundBestColumn);
    317             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    318           }
    319           else {
    320             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    321             solver->setColUpper(bestColumn, originalBoundBestColumn);
    322           }
    323         }
    324         else
    325           break;
    326       }
    327       else
    328         break;
    329     }
    330 
    331     if(!solver->isProvenOptimal())
    332       break;
    333 
    334     if(iteration > maxIterations_) {
    335       break;
    336     }
    337 
    338     if(CoinCpuTime()-time1 > maxTime_) {
    339       break;
    340     }
    341 
    342     memcpy(newSolution,solution,numberColumns*sizeof(double));
    343     numberFractionalVariables = 0;
    344     for (int i=0; i<numberIntegers; i++) {
    345       int iColumn = integerVariable[i];
    346       double value=newSolution[iColumn];
    347       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    348         numberFractionalVariables++;
    349       }
    350     }
    351 
    352   }
    353 
    354 
    355   double * rowActivity = new double[numberRows];
    356   memset(rowActivity,0,numberRows*sizeof(double));
    357 
    358   // re-compute new solution value
    359   double objOffset=0.0;
    360   solver->getDblParam(OsiObjOffset,objOffset);
    361   newSolutionValue = -objOffset;
    362   for (int i=0 ; i<numberColumns ; i++ )
    363     newSolutionValue += objective[i]*newSolution[i];
    364   newSolutionValue *= direction;
    365     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    366   if (newSolutionValue<solutionValue) {
    367     // paranoid check
    368     memset(rowActivity,0,numberRows*sizeof(double));
    369     for (int i=0;i<numberColumns;i++) {
    370       int j;
    371       double value = newSolution[i];
    372       if (value) {
    373         for (j=columnStart[i];
    374              j<columnStart[i]+columnLength[i];j++) {
    375           int iRow=row[j];
    376           rowActivity[iRow] += value*element[j];
    377         }
    378       }
    379     }
    380     // check was approximately feasible
    381     bool feasible=true;
    382     for (int i=0;i<numberRows;i++) {
    383       if(rowActivity[i]<rowLower[i]) {
    384         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    385           feasible = false;
    386       } else if(rowActivity[i]>rowUpper[i]) {
    387         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    388           feasible = false;
    389       }
    390     }
    391     for (int i=0; i<numberIntegers; i++) {
    392       int iColumn = integerVariable[i];
    393       double value=newSolution[iColumn];
    394       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    395         feasible = false;
    396         break;
    397       }
    398     }
    399     if (feasible) {
    400       // new solution
    401       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    402       solutionValue = newSolutionValue;
    403       //printf("** Solution of %g found by CbcHeuristicDiveFractional\n",newSolutionValue);
    404       returnCode=1;
    405     } else {
    406       // Can easily happen
    407       //printf("Debug CbcHeuristicDiveFractional giving bad solution\n");
    408     }
    409   }
    410 
    411   delete [] newSolution;
    412   delete [] columnFixed;
    413   delete [] originalBound;
    414   delete [] fixedAtLowerBound;
    415   delete [] rowActivity;
    416   delete solver;
    417   return returnCode;
    418 }
    419 
    420 // update model
    421 void CbcHeuristicDiveFractional::setModel(CbcModel * model)
    422 {
    423   model_ = model;
    424   // Get a copy of original matrix (and by row for rounding);
    425   assert(model_->solver());
    426   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    427   if (matrix) {
    428     matrix_ = *matrix;
    429     // make sure model okay for heuristic
    430     validate();
    431102  }
    432103}
    433 
    434 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    435 void
    436 CbcHeuristicDiveFractional::validate()
    437 {
    438   if (model_&&when()<10) {
    439     if (model_->numberIntegers()!=
    440         model_->numberObjects())
    441       setWhen(0);
    442   }
    443 
    444   int numberIntegers = model_->numberIntegers();
    445   const int * integerVariable = model_->integerVariable();
    446   delete [] downLocks_;
    447   delete [] upLocks_;
    448   downLocks_ = new unsigned short [numberIntegers];
    449   upLocks_ = new unsigned short [numberIntegers];
    450   // Column copy
    451   const double * element = matrix_.getElements();
    452   const int * row = matrix_.getIndices();
    453   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    454   const int * columnLength = matrix_.getVectorLengths();
    455   const double * rowLower = model_->solver()->getRowLower();
    456   const double * rowUpper = model_->solver()->getRowUpper();
    457   for (int i=0;i<numberIntegers;i++) {
    458     int iColumn = integerVariable[i];
    459     int down=0;
    460     int up=0;
    461     if (columnLength[iColumn]>65535) {
    462       setWhen(0);
    463       break; // unlikely to work
    464     }
    465     for (CoinBigIndex j=columnStart[iColumn];
    466          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    467       int iRow=row[j];
    468       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    469         up++;
    470         down++;
    471       } else if (element[j]>0.0) {
    472         if (rowUpper[iRow]<1.0e20)
    473           up++;
    474         else
    475           down++;
    476       } else {
    477         if (rowLower[iRow]>-1.0e20)
    478           up++;
    479         else
    480           down++;
    481       }
    482     }
    483     downLocks_[i] = (unsigned short) down;
    484     upLocks_[i] = (unsigned short) up;
    485   }
    486 }
Note: See TracChangeset for help on using the changeset viewer.