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

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