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

    r873 r912  
    66#endif
    77
     8//#define PRINT_DEBUG
     9
    810#include "CbcHeuristicDiveCoefficient.hpp"
    911#include "CbcStrategy.hpp"
    10 #include  "CoinTime.hpp"
    1112
    1213// Default Constructor
    1314CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient()
    14   :CbcHeuristic()
     15  :CbcHeuristicDive()
    1516{
    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;
    2217}
    2318
    2419// Constructor from model
    2520CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(CbcModel & model)
    26   :CbcHeuristic(model)
     21  :CbcHeuristicDive(model)
    2722{
    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;
    4023}
    4124
     
    4326CbcHeuristicDiveCoefficient::~CbcHeuristicDiveCoefficient ()
    4427{
    45   delete [] downLocks_;
    46   delete [] upLocks_;
    4728}
    4829
     
    6243  fprintf(fp,"3  CbcHeuristicDiveCoefficient heuristicDiveCoefficient(*cbcModel);\n");
    6344  CbcHeuristic::generateCpp(fp,"heuristicDiveCoefficient");
    64   if (percentageToFix_!=other.percentageToFix_)
    65     fprintf(fp,"3  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
    66   else
    67     fprintf(fp,"4  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
    68   if (maxIterations_!=other.maxIterations_)
    69     fprintf(fp,"3  heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);
    70   else
    71     fprintf(fp,"4  heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);
    72   if (maxTime_!=other.maxTime_)
    73     fprintf(fp,"3  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
    74   else
    75     fprintf(fp,"4  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
    7645  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveCoefficient);\n");
    7746}
     
    8049CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(const CbcHeuristicDiveCoefficient & rhs)
    8150:
    82   CbcHeuristic(rhs),
    83   matrix_(rhs.matrix_),
    84   percentageToFix_(rhs.percentageToFix_),
    85   maxIterations_(rhs.maxIterations_),
    86   maxTime_(rhs.maxTime_)
     51  CbcHeuristicDive(rhs)
    8752{
    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   }
    9753}
    9854
     
    10258{
    10359  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     }
     60    CbcHeuristicDive::operator=(rhs);
    11961  }
    12062  return *this;
    12163}
    12264
    123 // Resets stuff if model changes
    124 void
    125 CbcHeuristicDiveCoefficient::resetModel(CbcModel * model)
     65void
     66CbcHeuristicDiveCoefficient::selectVariableToBranch(OsiSolverInterface* solver,
     67                                                    const double* newSolution,
     68                                                    int& bestColumn,
     69                                                    int& bestRound)
    12670{
    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 CbcHeuristicDiveCoefficient::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());
    16671  int numberIntegers = model_->numberIntegers();
    16772  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();
     73  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    17674
    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;
     75  bestColumn = -1;
     76  bestRound = -1; // -1 rounds down, +1 rounds up
     77  double bestFraction = DBL_MAX;
     78  int bestLocks = COIN_INT_MAX;
    19179  for (int i=0; i<numberIntegers; i++) {
    19280    int iColumn = integerVariable[i];
    19381    double value=newSolution[iColumn];
     82    double fraction=value-floor(value);
     83    int round = 0;
    19484    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 bestLocks = COIN_INT_MAX;
    208     int numberAtBoundFixed = 0;
    209     bool canRoundSolution = true;
    210     for (int i=0; i<numberIntegers; i++) {
    211       int iColumn = integerVariable[i];
    212       double value=newSolution[iColumn];
    213       double fraction=value-floor(value);
    214       int round = 0;
    215       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    216         int nDownLocks = downLocks_[i];
    217         int nUpLocks = upLocks_[i];
    218         if(nDownLocks>0 && nUpLocks>0) {
    219           // the variable cannot be rounded
    220           canRoundSolution = false;
    221           int nLocks = nDownLocks;
    222           if(nDownLocks<nUpLocks)
    223             round = -1;
    224           else if(nDownLocks>nUpLocks) {
    225             round = 1;
    226             fraction = 1.0 - fraction;
    227             nLocks = nUpLocks;
    228           }
    229           else if(fraction < 0.5)
    230             round = -1;
    231           else {
    232             round = 1;
    233             fraction = 1.0 - fraction;
    234             nLocks = nUpLocks;
    235           }
    236 
    237           // if variable is not binary, penalize it
    238           if(!solver->isBinary(iColumn))
    239             fraction *= 1000.0;
    240 
    241           if(nLocks < bestLocks || (nLocks == bestLocks &&
    242                                     fraction < bestFraction)) {
    243             bestColumn = iColumn;
    244             bestLocks = nLocks;
    245             bestFraction = fraction;
    246             bestRound = round;
    247           }
     85      int nDownLocks = downLocks_[i];
     86      int nUpLocks = upLocks_[i];
     87      if(nDownLocks>0 && nUpLocks>0) {
     88        // the variable cannot be rounded
     89        int nLocks = nDownLocks;
     90        if(nDownLocks<nUpLocks)
     91          round = -1;
     92        else if(nDownLocks>nUpLocks) {
     93          round = 1;
     94          fraction = 1.0 - fraction;
     95          nLocks = nUpLocks;
    24896        }
    249       }
    250       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    251         // fix the variable at one of its bounds
    252         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    253             lower[iColumn] != upper[iColumn]) {
    254           columnFixed[numberAtBoundFixed] = iColumn;
    255           originalBound[numberAtBoundFixed] = upper[iColumn];
    256           fixedAtLowerBound[numberAtBoundFixed] = true;
    257           solver->setColUpper(iColumn, lower[iColumn]);
    258           numberAtBoundFixed++;
     97        else if(fraction < 0.5)
     98          round = -1;
     99        else {
     100          round = 1;
     101          fraction = 1.0 - fraction;
     102          nLocks = nUpLocks;
    259103        }
    260         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    261             lower[iColumn] != upper[iColumn]) {
    262           columnFixed[numberAtBoundFixed] = iColumn;
    263           originalBound[numberAtBoundFixed] = lower[iColumn];
    264           fixedAtLowerBound[numberAtBoundFixed] = false;
    265           solver->setColLower(iColumn, upper[iColumn]);
    266           numberAtBoundFixed++;
     104       
     105        // if variable is not binary, penalize it
     106        if(!solver->isBinary(iColumn))
     107          fraction *= 1000.0;
     108       
     109        if(nLocks < bestLocks || (nLocks == bestLocks &&
     110                                  fraction < bestFraction)) {
     111          bestColumn = iColumn;
     112          bestLocks = nLocks;
     113          bestFraction = fraction;
     114          bestRound = round;
    267115        }
    268116      }
    269117    }
    270 
    271     if(canRoundSolution) {
    272       // Round all the fractional variables
    273       for (int i=0; i<numberIntegers; i++) {
    274         int iColumn = integerVariable[i];
    275         double value=newSolution[iColumn];
    276         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    277           if(downLocks_[i]==0 || upLocks_[i]==0) {
    278             if(downLocks_[i]==0 && upLocks_[i]==0) {
    279               if(direction * objective[iColumn] >= 0.0)
    280                 newSolution[iColumn] = floor(value);
    281               else
    282                 newSolution[iColumn] = ceil(value);
    283             }
    284             else if(downLocks_[i]==0)
    285               newSolution[iColumn] = floor(value);
    286             else
    287               newSolution[iColumn] = ceil(value);
    288           }
    289           else
    290             break;
    291         }
    292       }
    293       break;
    294     }
    295 
    296 
    297     double originalBoundBestColumn;
    298     if(bestColumn >= 0) {
    299       if(bestRound < 0) {
    300         originalBoundBestColumn = upper[bestColumn];
    301         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    302       }
    303       else {
    304         originalBoundBestColumn = lower[bestColumn];
    305         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    306       }
    307     } else {
    308       break;
    309     }
    310     int originalBestRound = bestRound;
    311     while (1) {
    312 
    313       solver->resolve();
    314 
    315       if(!solver->isProvenOptimal()) {
    316         if(numberAtBoundFixed > 0) {
    317           // Remove the bound fix for variables that were at bounds
    318           for(int i=0; i<numberAtBoundFixed; i++) {
    319             int iColFixed = columnFixed[i];
    320             if(fixedAtLowerBound[i])
    321               solver->setColUpper(iColFixed, originalBound[i]);
    322             else
    323               solver->setColLower(iColFixed, originalBound[i]);
    324           }
    325           numberAtBoundFixed = 0;
    326         }
    327         else if(bestRound == originalBestRound) {
    328           bestRound *= (-1);
    329           if(bestRound < 0) {
    330             solver->setColLower(bestColumn, originalBoundBestColumn);
    331             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    332           }
    333           else {
    334             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    335             solver->setColUpper(bestColumn, originalBoundBestColumn);
    336           }
    337         }
    338         else
    339           break;
    340       }
    341       else
    342         break;
    343     }
    344 
    345     if(!solver->isProvenOptimal())
    346       break;
    347 
    348     if(iteration > maxIterations_) {
    349       break;
    350     }
    351 
    352     if(CoinCpuTime()-time1 > maxTime_) {
    353       break;
    354     }
    355 
    356     memcpy(newSolution,solution,numberColumns*sizeof(double));
    357     numberFractionalVariables = 0;
    358     for (int i=0; i<numberIntegers; i++) {
    359       int iColumn = integerVariable[i];
    360       double value=newSolution[iColumn];
    361       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    362         numberFractionalVariables++;
    363       }
    364     }
    365 
    366   }
    367 
    368 
    369   double * rowActivity = new double[numberRows];
    370   memset(rowActivity,0,numberRows*sizeof(double));
    371 
    372   // re-compute new solution value
    373   double objOffset=0.0;
    374   solver->getDblParam(OsiObjOffset,objOffset);
    375   newSolutionValue = -objOffset;
    376   for (int i=0 ; i<numberColumns ; i++ )
    377     newSolutionValue += objective[i]*newSolution[i];
    378   newSolutionValue *= direction;
    379     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    380   if (newSolutionValue<solutionValue) {
    381     // paranoid check
    382     memset(rowActivity,0,numberRows*sizeof(double));
    383     for (int i=0;i<numberColumns;i++) {
    384       int j;
    385       double value = newSolution[i];
    386       if (value) {
    387         for (j=columnStart[i];
    388              j<columnStart[i]+columnLength[i];j++) {
    389           int iRow=row[j];
    390           rowActivity[iRow] += value*element[j];
    391         }
    392       }
    393     }
    394     // check was approximately feasible
    395     bool feasible=true;
    396     for (int i=0;i<numberRows;i++) {
    397       if(rowActivity[i]<rowLower[i]) {
    398         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    399           feasible = false;
    400       } else if(rowActivity[i]>rowUpper[i]) {
    401         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    402           feasible = false;
    403       }
    404     }
    405     for (int i=0; i<numberIntegers; i++) {
    406       int iColumn = integerVariable[i];
    407       double value=newSolution[iColumn];
    408       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    409         feasible = false;
    410         break;
    411       }
    412     }
    413     if (feasible) {
    414       // new solution
    415       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    416       solutionValue = newSolutionValue;
    417       //printf("** Solution of %g found by CbcHeuristicDiveCoefficient\n",newSolutionValue);
    418       returnCode=1;
    419     } else {
    420       // Can easily happen
    421       //printf("Debug CbcHeuristicDiveCoefficient giving bad solution\n");
    422     }
    423   }
    424 
    425   delete [] newSolution;
    426   delete [] columnFixed;
    427   delete [] originalBound;
    428   delete [] fixedAtLowerBound;
    429   delete [] rowActivity;
    430   delete solver;
    431   return returnCode;
    432 }
    433 
    434 // update model
    435 void CbcHeuristicDiveCoefficient::setModel(CbcModel * model)
    436 {
    437   model_ = model;
    438   // Get a copy of original matrix (and by row for rounding);
    439   assert(model_->solver());
    440   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    441   if (matrix) {
    442     matrix_ = *matrix;
    443     // make sure model okay for heuristic
    444     validate();
    445118  }
    446119}
    447 
    448 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    449 void
    450 CbcHeuristicDiveCoefficient::validate()
    451 {
    452   if (model_&&when()<10) {
    453     if (model_->numberIntegers()!=
    454         model_->numberObjects())
    455       setWhen(0);
    456   }
    457 
    458   int numberIntegers = model_->numberIntegers();
    459   const int * integerVariable = model_->integerVariable();
    460   delete [] downLocks_;
    461   delete [] upLocks_;
    462   downLocks_ = new unsigned short [numberIntegers];
    463   upLocks_ = new unsigned short [numberIntegers];
    464   // Column copy
    465   const double * element = matrix_.getElements();
    466   const int * row = matrix_.getIndices();
    467   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    468   const int * columnLength = matrix_.getVectorLengths();
    469   const double * rowLower = model_->solver()->getRowLower();
    470   const double * rowUpper = model_->solver()->getRowUpper();
    471   for (int i=0;i<numberIntegers;i++) {
    472     int iColumn = integerVariable[i];
    473     int down=0;
    474     int up=0;
    475     if (columnLength[iColumn]>65535) {
    476       setWhen(0);
    477       break; // unlikely to work
    478     }
    479     for (CoinBigIndex j=columnStart[iColumn];
    480          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    481       int iRow=row[j];
    482       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    483         up++;
    484         down++;
    485       } else if (element[j]>0.0) {
    486         if (rowUpper[iRow]<1.0e20)
    487           up++;
    488         else
    489           down++;
    490       } else {
    491         if (rowLower[iRow]>-1.0e20)
    492           up++;
    493         else
    494           down++;
    495       }
    496     }
    497     downLocks_[i] = (unsigned short) down;
    498     upLocks_[i] = (unsigned short) up;
    499   }
    500 }
Note: See TracChangeset for help on using the changeset viewer.