Changeset 906


Ignore:
Timestamp:
Apr 9, 2008 3:02:13 PM (11 years ago)
Author:
jpgoncal
Message:

Merged with changes done in trunk

Location:
branches/heur/Cbc/src
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • branches/heur/Cbc/src/CbcHeuristic.cpp

    r901 r906  
    1212#include <cfloat>
    1313
    14 #define PRINT_DEBUG
     14//#define PRINT_DEBUG
    1515#ifdef COIN_HAS_CLP
    1616#include "OsiClpSolverInterface.hpp"
     
    228228}
    229229
     230void
     231CbcHeuristic::printDistanceToNodes()
     232{
     233  const CbcNode* currentNode = model_->currentNode();
     234  if (currentNode != NULL) {
     235    CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
     236    for (int i = runNodes_.size() - 1; i >= 0; --i) {
     237      nodeDesc->distance(runNodes_.node(i));
     238    }
     239    runNodes_.append(nodeDesc);
     240  }
     241}
     242
    230243bool
    231244CbcHeuristic::shouldHeurRun()
     
    295308    // Get where we are and create the appropriate CbcHeuristicNode object
    296309    CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
    297 #ifdef PRINT_DEBUG
     310    //#ifdef PRINT_DEBUG
     311#if 1
    298312    double minDistance = nodeDesc->minDistance(runNodes_);
    299     double minDistanceToRun = log(depth) / log(2);
     313    double minDistanceToRun = 1.5 * log(depth) / log(2);
    300314    std::cout<<"minDistance = "<<minDistance
    301315             <<", minDistanceToRun = "<<minDistanceToRun<<std::endl;
     
    326340  if(depth != 0) {
    327341
    328     double probability = depth / pow(2,depth);
     342    double probability = pow(depth,2) / pow(2,depth);
    329343    double randomNumber = randomNumberGenerator_.randomDouble();
    330344    if (randomNumber>probability)
     
    812826  const double overlapWeight = 0.4;
    813827  const double subsetWeight = 0.2;
     828  int countDisjointWeight = 0;
     829  int countOverlapWeight = 0;
     830  int countSubsetWeight = 0;
    814831  int i = 0;
    815832  int j = 0;
     
    845862    if (brComp < 0) {
    846863      dist += subsetWeight;
     864      countSubsetWeight++;
    847865      ++i;
    848866    }
    849867    else if (brComp > 0) {
    850868      dist += subsetWeight;
     869      countSubsetWeight++;
    851870      ++j;
    852871    }
     
    859878      case CbcRangeDisjoint: // disjoint decisions
    860879        dist += disjointWeight;
     880        countDisjointWeight++;
    861881        break;
    862882      case CbcRangeSubset: // subset one way or another
    863883      case CbcRangeSuperset:
    864884        dist += subsetWeight;
     885        countSubsetWeight++;
    865886        break;
    866887      case CbcRangeOverlap: // overlap
    867888        dist += overlapWeight;
     889        countOverlapWeight++;
    868890        break;
    869891      }
     
    873895  }
    874896  dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
     897  countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
     898  printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
     899         countOverlapWeight, countDisjointWeight);
    875900  return dist;
    876901}
  • branches/heur/Cbc/src/CbcHeuristic.hpp

    r901 r906  
    181181  bool shouldHeurRun_randomChoice();
    182182  void debugNodes();
     183  void printDistanceToNodes();
    183184
    184185protected:
  • branches/heur/Cbc/src/CbcHeuristicDiveCoefficient.cpp

    r901 r906  
    66#endif
    77
    8 #define PRINT_DEBUG
     8//#define PRINT_DEBUG
    99
    1010#include "CbcHeuristicDiveCoefficient.hpp"
    1111#include "CbcStrategy.hpp"
    12 #include  "CoinTime.hpp"
    1312
    1413// Default Constructor
    1514CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient()
    16   :CbcHeuristic()
     15  :CbcHeuristicDive()
    1716{
    18   // matrix and row copy will automatically be empty
    19   downLocks_ =NULL;
    20   upLocks_ = NULL;
    21   percentageToFix_ = 0.2;
    22   maxIterations_ = 100;
    23   maxTime_ = 60;
    2417}
    2518
    2619// Constructor from model
    2720CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(CbcModel & model)
    28   :CbcHeuristic(model)
     21  :CbcHeuristicDive(model)
    2922{
    30   downLocks_ =NULL;
    31   upLocks_ = NULL;
    32   // Get a copy of original matrix
    33   assert(model.solver());
    34   // model may have empty matrix - wait until setModel
    35   const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
    36   if (matrix) {
    37     matrix_ = *matrix;
    38     validate();
    39   }
    40   percentageToFix_ = 0.2;
    41   maxIterations_ = 100;
    42   maxTime_ = 60;
    4323}
    4424
     
    4626CbcHeuristicDiveCoefficient::~CbcHeuristicDiveCoefficient ()
    4727{
    48   delete [] downLocks_;
    49   delete [] upLocks_;
    5028}
    5129
     
    6543  fprintf(fp,"3  CbcHeuristicDiveCoefficient heuristicDiveCoefficient(*cbcModel);\n");
    6644  CbcHeuristic::generateCpp(fp,"heuristicDiveCoefficient");
    67   if (percentageToFix_!=other.percentageToFix_)
    68     fprintf(fp,"3  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
    69   else
    70     fprintf(fp,"4  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
    71   if (maxIterations_!=other.maxIterations_)
    72     fprintf(fp,"3  heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);
    73   else
    74     fprintf(fp,"4  heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);
    75   if (maxTime_!=other.maxTime_)
    76     fprintf(fp,"3  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
    77   else
    78     fprintf(fp,"4  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
    7945  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveCoefficient);\n");
    8046}
     
    8349CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(const CbcHeuristicDiveCoefficient & rhs)
    8450:
    85   CbcHeuristic(rhs),
    86   matrix_(rhs.matrix_),
    87   percentageToFix_(rhs.percentageToFix_),
    88   maxIterations_(rhs.maxIterations_),
    89   maxTime_(rhs.maxTime_)
     51  CbcHeuristicDive(rhs)
    9052{
    91   int numberIntegers = model_->numberIntegers();
    92   if (rhs.downLocks_) {
    93     int numberIntegers = model_->numberIntegers();
    94     downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    95     upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    96   } else {
    97     downLocks_ = NULL;
    98     upLocks_ = NULL;
    99   }
    10053}
    10154
     
    10558{
    10659  if (this!=&rhs) {
    107     CbcHeuristic::operator=(rhs);
    108     matrix_ = rhs.matrix_;
    109     percentageToFix_ = rhs.percentageToFix_;
    110     maxIterations_ = rhs.maxIterations_;
    111     maxTime_ = rhs.maxTime_;
    112     delete [] downLocks_;
    113     delete [] upLocks_;
    114     if (rhs.downLocks_) {
    115       int numberIntegers = model_->numberIntegers();
    116       downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
    117       upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
    118     } else {
    119       downLocks_ = NULL;
    120       upLocks_ = NULL;
    121     }
     60    CbcHeuristicDive::operator=(rhs);
    12261  }
    12362  return *this;
    12463}
    12564
    126 // Resets stuff if model changes
    127 void
    128 CbcHeuristicDiveCoefficient::resetModel(CbcModel * model)
     65void
     66CbcHeuristicDiveCoefficient::selectVariableToBranch(OsiSolverInterface* solver,
     67                                                    const double* newSolution,
     68                                                    int& bestColumn,
     69                                                    int& bestRound)
    12970{
    130   model_=model;
    131   // Get a copy of original matrix
    132   assert(model_->solver());
    133   // model may have empty matrix - wait until setModel
    134   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    135   if (matrix) {
    136     matrix_ = *matrix;
    137     validate();
    138   }
    139 }
    140 
    141 // See if dive fractional will give better solution
    142 // Sets value of solution
    143 // Returns 1 if solution, 0 if not
    144 int
    145 CbcHeuristicDiveCoefficient::solution(double & solutionValue,
    146                                       double * betterSolution)
    147 {
    148   ++numCouldRun_;
    149 
    150   //  bool runHeuristic = shouldHeurRun();
    151   bool runHeuristic = shouldHeurRun_randomChoice();
    152   //  std::cout<<"shouldHeurRun = "<<runHeuristic<<std::endl;
    153 
    154 #ifdef PRINT_DEBUG
    155   std::cout<<"numCouldRun_ = "<<numCouldRun_
    156            <<", numRuns_ = "<<numRuns_;
    157 #endif
    158 
    159 
    160   if (! runHeuristic) {
    161 #ifdef PRINT_DEBUG
    162     std::cout<<", solutionFound = 0"<<std::endl;
    163 #endif
    164     return 0;
    165   }
    166 
    167 #if 0
    168   // See if to do
    169   if (!when()||(when()%10==1&&model_->phase()!=1)||
    170       (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
    171     return 0; // switched off
    172 #endif
    173 
    174   double time1 = CoinCpuTime();
    175 
    176   OsiSolverInterface * solver = model_->solver()->clone();
    177   const double * lower = solver->getColLower();
    178   const double * upper = solver->getColUpper();
    179   const double * rowLower = solver->getRowLower();
    180   const double * rowUpper = solver->getRowUpper();
    181   const double * solution = solver->getColSolution();
    182   const double * objective = solver->getObjCoefficients();
    183   double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    184   double primalTolerance;
    185   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    186 
    187   int numberRows = matrix_.getNumRows();
    188   assert (numberRows<=solver->getNumRows());
    18971  int numberIntegers = model_->numberIntegers();
    19072  const int * integerVariable = model_->integerVariable();
    191   double direction = solver->getObjSense(); // 1 for min, -1 for max
    192   double newSolutionValue = direction*solver->getObjValue();
    193   int returnCode = 0;
    194   // Column copy
    195   const double * element = matrix_.getElements();
    196   const int * row = matrix_.getIndices();
    197   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    198   const int * columnLength = matrix_.getVectorLengths();
     73  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    19974
    200   // Get solution array for heuristic solution
    201   int numberColumns = solver->getNumCols();
    202   double * newSolution = new double [numberColumns];
    203   memcpy(newSolution,solution,numberColumns*sizeof(double));
    204 
    205   // vectors to store the latest variables fixed at their bounds
    206   int* columnFixed = new int [numberIntegers];
    207   double* originalBound = new double [numberIntegers];
    208   bool * fixedAtLowerBound = new bool [numberIntegers];
    209 
    210   const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
    211 
    212   // count how many fractional variables
    213   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;
    21479  for (int i=0; i<numberIntegers; i++) {
    21580    int iColumn = integerVariable[i];
    21681    double value=newSolution[iColumn];
     82    double fraction=value-floor(value);
     83    int round = 0;
    21784    if (fabs(floor(value+0.5)-value)>integerTolerance) {
    218       numberFractionalVariables++;
    219     }
    220   }
    221 
    222   int iteration = 0;
    223   while(numberFractionalVariables) {
    224     iteration++;
    225 
    226     // select a fractional variable to bound
    227     double bestFraction = DBL_MAX;
    228     int bestColumn = -1;
    229     int bestRound = -1; // -1 rounds down, +1 rounds up
    230     int bestLocks = COIN_INT_MAX;
    231     int numberAtBoundFixed = 0;
    232     bool canRoundSolution = true;
    233     for (int i=0; i<numberIntegers; i++) {
    234       int iColumn = integerVariable[i];
    235       double value=newSolution[iColumn];
    236       double fraction=value-floor(value);
    237       int round = 0;
    238       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    239         int nDownLocks = downLocks_[i];
    240         int nUpLocks = upLocks_[i];
    241         if(nDownLocks>0 && nUpLocks>0) {
    242           // the variable cannot be rounded
    243           canRoundSolution = false;
    244           int nLocks = nDownLocks;
    245           if(nDownLocks<nUpLocks)
    246             round = -1;
    247           else if(nDownLocks>nUpLocks) {
    248             round = 1;
    249             fraction = 1.0 - fraction;
    250             nLocks = nUpLocks;
    251           }
    252           else if(fraction < 0.5)
    253             round = -1;
    254           else {
    255             round = 1;
    256             fraction = 1.0 - fraction;
    257             nLocks = nUpLocks;
    258           }
    259 
    260           // if variable is not binary, penalize it
    261           if(!solver->isBinary(iColumn))
    262             fraction *= 1000.0;
    263 
    264           if(nLocks < bestLocks || (nLocks == bestLocks &&
    265                                     fraction < bestFraction)) {
    266             bestColumn = iColumn;
    267             bestLocks = nLocks;
    268             bestFraction = fraction;
    269             bestRound = round;
    270           }
     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;
    27196        }
    272       }
    273       else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
    274         // fix the variable at one of its bounds
    275         if (fabs(lower[iColumn]-value)<=integerTolerance &&
    276             lower[iColumn] != upper[iColumn]) {
    277           columnFixed[numberAtBoundFixed] = iColumn;
    278           originalBound[numberAtBoundFixed] = upper[iColumn];
    279           fixedAtLowerBound[numberAtBoundFixed] = true;
    280           solver->setColUpper(iColumn, lower[iColumn]);
    281           numberAtBoundFixed++;
     97        else if(fraction < 0.5)
     98          round = -1;
     99        else {
     100          round = 1;
     101          fraction = 1.0 - fraction;
     102          nLocks = nUpLocks;
    282103        }
    283         else if(fabs(upper[iColumn]-value)<=integerTolerance &&
    284             lower[iColumn] != upper[iColumn]) {
    285           columnFixed[numberAtBoundFixed] = iColumn;
    286           originalBound[numberAtBoundFixed] = lower[iColumn];
    287           fixedAtLowerBound[numberAtBoundFixed] = false;
    288           solver->setColLower(iColumn, upper[iColumn]);
    289           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;
    290115        }
    291116      }
    292117    }
    293 
    294     if(canRoundSolution) {
    295       // Round all the fractional variables
    296       for (int i=0; i<numberIntegers; i++) {
    297         int iColumn = integerVariable[i];
    298         double value=newSolution[iColumn];
    299         if (fabs(floor(value+0.5)-value)>integerTolerance) {
    300           if(downLocks_[i]==0 || upLocks_[i]==0) {
    301             if(downLocks_[i]==0 && upLocks_[i]==0) {
    302               if(direction * objective[iColumn] >= 0.0)
    303                 newSolution[iColumn] = floor(value);
    304               else
    305                 newSolution[iColumn] = ceil(value);
    306             }
    307             else if(downLocks_[i]==0)
    308               newSolution[iColumn] = floor(value);
    309             else
    310               newSolution[iColumn] = ceil(value);
    311           }
    312           else
    313             break;
    314         }
    315       }
    316       break;
    317     }
    318 
    319 
    320     double originalBoundBestColumn;
    321     if(bestColumn >= 0) {
    322       if(bestRound < 0) {
    323         originalBoundBestColumn = upper[bestColumn];
    324         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    325       }
    326       else {
    327         originalBoundBestColumn = lower[bestColumn];
    328         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    329       }
    330     } else {
    331       break;
    332     }
    333     int originalBestRound = bestRound;
    334     while (1) {
    335 
    336       solver->resolve();
    337 
    338       if(!solver->isProvenOptimal()) {
    339         if(numberAtBoundFixed > 0) {
    340           // Remove the bound fix for variables that were at bounds
    341           for(int i=0; i<numberAtBoundFixed; i++) {
    342             int iColFixed = columnFixed[i];
    343             if(fixedAtLowerBound[i])
    344               solver->setColUpper(iColFixed, originalBound[i]);
    345             else
    346               solver->setColLower(iColFixed, originalBound[i]);
    347           }
    348           numberAtBoundFixed = 0;
    349         }
    350         else if(bestRound == originalBestRound) {
    351           bestRound *= (-1);
    352           if(bestRound < 0) {
    353             solver->setColLower(bestColumn, originalBoundBestColumn);
    354             solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
    355           }
    356           else {
    357             solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
    358             solver->setColUpper(bestColumn, originalBoundBestColumn);
    359           }
    360         }
    361         else
    362           break;
    363       }
    364       else
    365         break;
    366     }
    367 
    368     if(!solver->isProvenOptimal())
    369       break;
    370 
    371     if(iteration > maxIterations_) {
    372       break;
    373     }
    374 
    375     if(CoinCpuTime()-time1 > maxTime_) {
    376       break;
    377     }
    378 
    379     memcpy(newSolution,solution,numberColumns*sizeof(double));
    380     numberFractionalVariables = 0;
    381     for (int i=0; i<numberIntegers; i++) {
    382       int iColumn = integerVariable[i];
    383       double value=newSolution[iColumn];
    384       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    385         numberFractionalVariables++;
    386       }
    387     }
    388 
    389   }
    390 
    391   bool updateHowOften = false;
    392   // Good question when... every time the heur was run? every time it
    393   // has found a solution? every time it has found a better solution?
    394   // for now update every time
    395   updateHowOften = true;
    396 
    397   double * rowActivity = new double[numberRows];
    398   memset(rowActivity,0,numberRows*sizeof(double));
    399 
    400   // re-compute new solution value
    401   double objOffset=0.0;
    402   solver->getDblParam(OsiObjOffset,objOffset);
    403   newSolutionValue = -objOffset;
    404   for (int i=0 ; i<numberColumns ; i++ )
    405     newSolutionValue += objective[i]*newSolution[i];
    406   newSolutionValue *= direction;
    407     //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    408   if (newSolutionValue<solutionValue) {
    409     // paranoid check
    410     memset(rowActivity,0,numberRows*sizeof(double));
    411     for (int i=0;i<numberColumns;i++) {
    412       int j;
    413       double value = newSolution[i];
    414       if (value) {
    415         for (j=columnStart[i];
    416              j<columnStart[i]+columnLength[i];j++) {
    417           int iRow=row[j];
    418           rowActivity[iRow] += value*element[j];
    419         }
    420       }
    421     }
    422     // check was approximately feasible
    423     bool feasible=true;
    424     for (int i=0;i<numberRows;i++) {
    425       if(rowActivity[i]<rowLower[i]) {
    426         if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
    427           feasible = false;
    428       } else if(rowActivity[i]>rowUpper[i]) {
    429         if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
    430           feasible = false;
    431       }
    432     }
    433     for (int i=0; i<numberIntegers; i++) {
    434       int iColumn = integerVariable[i];
    435       double value=newSolution[iColumn];
    436       if (fabs(floor(value+0.5)-value)>integerTolerance) {
    437         feasible = false;
    438         break;
    439       }
    440     }
    441     if (feasible) {
    442       // new solution
    443       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    444       solutionValue = newSolutionValue;
    445       //printf("** Solution of %g found by CbcHeuristicDiveCoefficient\n",newSolutionValue);
    446       returnCode=1;
    447 
    448       //      setCurrentNode(model_->solver());
    449     }
    450   }
    451 
    452 #if 0
    453   if (updateHowOften) {
    454     howOften_ += (int) (howOften_*decayFactor_);
    455   }
    456 #endif
    457 
    458 #ifdef PRINT_DEBUG
    459   std::cout<<", solutionFound = "<<returnCode<<std::endl;
    460 #endif
    461 
    462   delete [] newSolution;
    463   delete [] columnFixed;
    464   delete [] originalBound;
    465   delete [] fixedAtLowerBound;
    466   delete [] rowActivity;
    467   delete solver;
    468   return returnCode;
    469 }
    470 
    471 // update model
    472 void CbcHeuristicDiveCoefficient::setModel(CbcModel * model)
    473 {
    474   model_ = model;
    475   // Get a copy of original matrix (and by row for rounding);
    476   assert(model_->solver());
    477   const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
    478   if (matrix) {
    479     matrix_ = *matrix;
    480     // make sure model okay for heuristic
    481     validate();
    482118  }
    483119}
    484 
    485 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    486 void
    487 CbcHeuristicDiveCoefficient::validate()
    488 {
    489   if (model_&&when()<10) {
    490     if (model_->numberIntegers()!=
    491         model_->numberObjects())
    492       setWhen(0);
    493   }
    494 
    495   int numberIntegers = model_->numberIntegers();
    496   const int * integerVariable = model_->integerVariable();
    497   delete [] downLocks_;
    498   delete [] upLocks_;
    499   downLocks_ = new unsigned short [numberIntegers];
    500   upLocks_ = new unsigned short [numberIntegers];
    501   // Column copy
    502   const double * element = matrix_.getElements();
    503   const int * row = matrix_.getIndices();
    504   const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    505   const int * columnLength = matrix_.getVectorLengths();
    506   const double * rowLower = model_->solver()->getRowLower();
    507   const double * rowUpper = model_->solver()->getRowUpper();
    508   for (int i=0;i<numberIntegers;i++) {
    509     int iColumn = integerVariable[i];
    510     int down=0;
    511     int up=0;
    512     if (columnLength[iColumn]>65535) {
    513       setWhen(0);
    514       break; // unlikely to work
    515     }
    516     for (CoinBigIndex j=columnStart[iColumn];
    517          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    518       int iRow=row[j];
    519       if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
    520         up++;
    521         down++;
    522       } else if (element[j]>0.0) {
    523         if (rowUpper[iRow]<1.0e20)
    524           up++;
    525         else
    526           down++;
    527       } else {
    528         if (rowLower[iRow]>-1.0e20)
    529           up++;
    530         else
    531           down++;
    532       }
    533     }
    534     downLocks_[i] = (unsigned short) down;
    535     upLocks_[i] = (unsigned short) up;
    536   }
    537 }
  • branches/heur/Cbc/src/CbcHeuristicDiveCoefficient.hpp

    r868 r906  
    44#define CbcHeuristicDiveCoefficient_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveCoefficient class
    99 */
    1010
    11 class CbcHeuristicDiveCoefficient : public CbcHeuristic {
     11class CbcHeuristicDiveCoefficient : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
    37 
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does Coefficient Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     35  /// Selects the next variable to branch on
     36  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     37                                      const double* newSolution,
     38                                      int& bestColumn,
     39                                      int& bestRound);
    8640
    8741};
  • branches/heur/Cbc/src/CbcHeuristicDiveFractional.cpp

    r871 r906  
    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 }
  • branches/heur/Cbc/src/CbcHeuristicDiveFractional.hpp

    r868 r906  
    44#define CbcHeuristicDiveFractional_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveFractional class
    99 */
    1010
    11 class CbcHeuristicDiveFractional : public CbcHeuristic {
     11class CbcHeuristicDiveFractional : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
    37 
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does Fractional Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     35  /// Selects the next variable to branch on
     36  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     37                                      const double* newSolution,
     38                                      int& bestColumn,
     39                                      int& bestRound);
    8640
    8741};
  • branches/heur/Cbc/src/CbcHeuristicDiveGuided.cpp

    r871 r906  
    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 }
  • branches/heur/Cbc/src/CbcHeuristicDiveGuided.hpp

    r868 r906  
    44#define CbcHeuristicDiveGuided_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveGuided class
    99 */
    1010
    11 class CbcHeuristicDiveGuided : public CbcHeuristic {
     11class CbcHeuristicDiveGuided : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
     35  /// Tests if the heuristic can run
     36  virtual bool canHeuristicRun();
    3737
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does Guided Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     38  /// Selects the next variable to branch on
     39  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     40                                      const double* newSolution,
     41                                      int& bestColumn,
     42                                      int& bestRound);
    8643
    8744};
  • branches/heur/Cbc/src/CbcHeuristicDiveVectorLength.cpp

    r871 r906  
    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 }
  • branches/heur/Cbc/src/CbcHeuristicDiveVectorLength.hpp

    r868 r906  
    44#define CbcHeuristicDiveVectorLength_H
    55
    6 #include "CbcHeuristic.hpp"
     6#include "CbcHeuristicDive.hpp"
    77
    88/** DiveVectorLength class
    99 */
    1010
    11 class CbcHeuristicDiveVectorLength : public CbcHeuristic {
     11class CbcHeuristicDiveVectorLength : public CbcHeuristicDive {
    1212public:
    1313
     
    3333  virtual void generateCpp( FILE * fp) ;
    3434
    35   /// Resets stuff if model changes
    36   virtual void resetModel(CbcModel * model);
    37 
    38   /// update model (This is needed if cliques update matrix etc)
    39   virtual void setModel(CbcModel * model);
    40  
    41   using CbcHeuristic::solution ;
    42   /** returns 0 if no solution, 1 if valid solution
    43       with better objective value than one passed in
    44       Sets solution values if good, sets objective value (only if good)
    45       This is called after cuts have been added - so can not add cuts
    46       This does VectorLength Diving
    47   */
    48   virtual int solution(double & objectiveValue,
    49                        double * newSolution);
    50 
    51   /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    52   virtual void validate();
    53 
    54   /// Set percentage of integer variables to fix at bounds
    55   void setPercentageToFix(double value)
    56   { percentageToFix_ = value; }
    57 
    58   /// Set maximum number of iterations
    59   void setMaxIterations(int value)
    60   { maxIterations_ = value; }
    61 
    62   /// Set maximum time allowed
    63   void setMaxTime(double value)
    64   { maxTime_ = value; }
    65 
    66 protected:
    67   // Data
    68 
    69   // Original matrix by column
    70   CoinPackedMatrix matrix_;
    71 
    72   // Down locks
    73   unsigned short * downLocks_;
    74 
    75   // Up locks
    76   unsigned short * upLocks_;
    77 
    78   // Percentage of integer variables to fix at bounds
    79   double percentageToFix_;
    80 
    81   // Maximum number of iterations
    82   int maxIterations_;
    83 
    84   // Maximum time allowed
    85   double maxTime_;
     35  /// Selects the next variable to branch on
     36  virtual void selectVariableToBranch(OsiSolverInterface* solver,
     37                                      const double* newSolution,
     38                                      int& bestColumn,
     39                                      int& bestRound);
    8640
    8741};
  • branches/heur/Cbc/src/Makefile.am

    r871 r906  
    3535        CbcFeasibilityBase.hpp \
    3636        CbcHeuristic.cpp CbcHeuristic.hpp \
     37        CbcHeuristicDive.cpp CbcHeuristicDive.hpp \
    3738        CbcHeuristicDiveCoefficient.cpp CbcHeuristicDiveCoefficient.hpp \
    3839        CbcHeuristicDiveFractional.cpp CbcHeuristicDiveFractional.hpp \
     
    307308        CbcFeasibilityBase.hpp \
    308309        CbcHeuristic.hpp \
     310        CbcHeuristicDive.hpp \
    309311        CbcHeuristicDiveCoefficient.hpp \
    310312        CbcHeuristicDiveFractional.hpp \
  • branches/heur/Cbc/src/Makefile.in

    r871 r906  
    162162        CbcCutGenerator.lo CbcEventHandler.lo CbcFathom.lo \
    163163        CbcFathomDynamicProgramming.lo CbcHeuristic.lo \
    164         CbcHeuristicDiveCoefficient.lo CbcHeuristicDiveFractional.lo \
    165         CbcHeuristicDiveGuided.lo CbcHeuristicDiveVectorLength.lo \
    166         CbcHeuristicFPump.lo CbcHeuristicGreedy.lo \
    167         CbcHeuristicLocal.lo CbcHeuristicRINS.lo CbcMessage.lo \
    168         CbcModel.lo CbcNode.lo CbcStatistics.lo CbcStrategy.lo \
    169         CbcTree.lo CbcTreeLocal.lo
     164        CbcHeuristicDive.lo CbcHeuristicDiveCoefficient.lo \
     165        CbcHeuristicDiveFractional.lo CbcHeuristicDiveGuided.lo \
     166        CbcHeuristicDiveVectorLength.lo CbcHeuristicFPump.lo \
     167        CbcHeuristicGreedy.lo CbcHeuristicLocal.lo CbcHeuristicRINS.lo \
     168        CbcMessage.lo CbcModel.lo CbcNode.lo CbcStatistics.lo \
     169        CbcStrategy.lo CbcTree.lo CbcTreeLocal.lo
    170170libCbc_la_OBJECTS = $(am_libCbc_la_OBJECTS)
    171171libCbcSolver_la_LIBADD =
     
    510510        CbcFeasibilityBase.hpp \
    511511        CbcHeuristic.cpp CbcHeuristic.hpp \
     512        CbcHeuristicDive.cpp CbcHeuristicDive.hpp \
    512513        CbcHeuristicDiveCoefficient.cpp CbcHeuristicDiveCoefficient.hpp \
    513514        CbcHeuristicDiveFractional.cpp CbcHeuristicDiveFractional.hpp \
     
    650651        CbcFeasibilityBase.hpp \
    651652        CbcHeuristic.hpp \
     653        CbcHeuristicDive.hpp \
    652654        CbcHeuristicDiveCoefficient.hpp \
    653655        CbcHeuristicDiveFractional.hpp \
     
    806808@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcGeneric.Po@am__quote@
    807809@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristic.Plo@am__quote@
     810@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristicDive.Plo@am__quote@
    808811@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristicDiveCoefficient.Plo@am__quote@
    809812@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcHeuristicDiveFractional.Plo@am__quote@
Note: See TracChangeset for help on using the changeset viewer.