Ignore:
Timestamp:
Nov 9, 2009 6:33:07 PM (10 years ago)
Author:
EdwinStraver
Message:

Changed formatting using AStyle -A4 -p

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/Cbc/src/CbcHeuristicPivotAndFix.cpp

    r1271 r1286  
    2222
    2323// Default Constructor
    24 CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix() 
    25   :CbcHeuristic()
     24CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix()
     25        : CbcHeuristic()
    2626{
    2727}
     
    3030
    3131CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix(CbcModel & model)
    32   :CbcHeuristic(model)
    33 {
    34 }
    35 
    36 // Destructor 
     32        : CbcHeuristic(model)
     33{
     34}
     35
     36// Destructor
    3737CbcHeuristicPivotAndFix::~CbcHeuristicPivotAndFix ()
    3838{
     
    4343CbcHeuristicPivotAndFix::clone() const
    4444{
    45   return new CbcHeuristicPivotAndFix(*this);
     45    return new CbcHeuristicPivotAndFix(*this);
    4646}
    4747// Create C++ lines to get to current state
    48 void 
    49 CbcHeuristicPivotAndFix::generateCpp( FILE * fp) 
    50 {
    51   CbcHeuristicPivotAndFix other;
    52   fprintf(fp,"0#include \"CbcHeuristicPivotAndFix.hpp\"\n");
    53   fprintf(fp,"3  CbcHeuristicPivotAndFix heuristicPFX(*cbcModel);\n");
    54   CbcHeuristic::generateCpp(fp,"heuristicPFX");
    55   fprintf(fp,"3  cbcModel->addHeuristic(&heuristicPFX);\n");
    56 }
    57 
    58 // Copy constructor 
     48void
     49CbcHeuristicPivotAndFix::generateCpp( FILE * fp)
     50{
     51    CbcHeuristicPivotAndFix other;
     52    fprintf(fp, "0#include \"CbcHeuristicPivotAndFix.hpp\"\n");
     53    fprintf(fp, "3  CbcHeuristicPivotAndFix heuristicPFX(*cbcModel);\n");
     54    CbcHeuristic::generateCpp(fp, "heuristicPFX");
     55    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicPFX);\n");
     56}
     57
     58// Copy constructor
    5959CbcHeuristicPivotAndFix::CbcHeuristicPivotAndFix(const CbcHeuristicPivotAndFix & rhs)
    60 :
    61   CbcHeuristic(rhs)
    62 {
    63 }
    64 
    65 // Assignment operator 
    66 CbcHeuristicPivotAndFix & 
    67 CbcHeuristicPivotAndFix::operator=( const CbcHeuristicPivotAndFix& rhs)
    68 {
    69   if (this!=&rhs) {
    70     CbcHeuristic::operator=(rhs);
    71   }
    72   return *this;
     60        :
     61        CbcHeuristic(rhs)
     62{
     63}
     64
     65// Assignment operator
     66CbcHeuristicPivotAndFix &
     67CbcHeuristicPivotAndFix::operator=( const CbcHeuristicPivotAndFix & rhs)
     68{
     69    if (this != &rhs) {
     70        CbcHeuristic::operator=(rhs);
     71    }
     72    return *this;
    7373}
    7474
    7575// Resets stuff if model changes
    76 void 
     76void
    7777CbcHeuristicPivotAndFix::resetModel(CbcModel * /*model*/)
    7878{
    79   //CbcHeuristic::resetModel(model);
     79    //CbcHeuristic::resetModel(model);
    8080}
    8181/*
     
    8484int
    8585CbcHeuristicPivotAndFix::solution(double & /*solutionValue*/,
    86                                   double * /*betterSolution*/)
    87 {
    88 
    89   numCouldRun_++; // Todo: Ask JJHF what this for.
    90   std::cout << "Entering Pivot-and-Fix Heuristic" << std::endl;
     86                                  double * /*betterSolution*/)
     87{
     88
     89    numCouldRun_++; // Todo: Ask JJHF what this for.
     90    std::cout << "Entering Pivot-and-Fix Heuristic" << std::endl;
    9191
    9292#ifdef FORNOW
    93   std::cout << "Lucky you! You're in the Pivot-and-Fix Heuristic" << std::endl;
    94   // The struct should be moved to member data
    95 
    96 
    97 
    98   typedef struct {
    99     int numberSolutions;
    100     int maximumSolutions;
    101     int numberColumns;
    102     double ** solution;
    103     int * numberUnsatisfied;
    104   } clpSolution;
    105 
    106   double start = CoinCpuTime();
    107 
    108   OsiClpSolverInterface * clpSolverOriginal
     93    std::cout << "Lucky you! You're in the Pivot-and-Fix Heuristic" << std::endl;
     94    // The struct should be moved to member data
     95
     96
     97
     98    typedef struct {
     99        int numberSolutions;
     100        int maximumSolutions;
     101        int numberColumns;
     102        double ** solution;
     103        int * numberUnsatisfied;
     104    } clpSolution;
     105
     106    double start = CoinCpuTime();
     107
     108    OsiClpSolverInterface * clpSolverOriginal
    109109    = dynamic_cast<OsiClpSolverInterface *> (model_->solver());
    110   assert (clpSolverOriginal);
    111   OsiClpSolverInterface *clpSolver(clpSolverOriginal);
    112  
    113   ClpSimplex * simplex = clpSolver->getModelPtr();
    114 
    115   // Initialize the structure holding the solutions
    116   clpSolution solutions;
    117   // Set typeStruct field of ClpTrustedData struct to one.
    118   // This tells Clp it's "Mahdi!"
    119   ClpTrustedData trustedSolutions;
    120   trustedSolutions.typeStruct = 1;
    121   trustedSolutions.data = &solutions;
    122   solutions.numberSolutions=0;
    123   solutions.maximumSolutions=0;
    124   solutions.numberColumns=simplex->numberColumns();
    125   solutions.solution=NULL;
    126   solutions.numberUnsatisfied=NULL;
    127   simplex->setTrustedUserPointer(&trustedSolutions);
    128                                          
    129   // Solve from all slack to get some points
    130   simplex->allSlackBasis();
    131   simplex->primal();
    132 
    133   // -------------------------------------------------
    134   // Get the problem information
    135   // - get the number of cols and rows
    136   int numCols = clpSolver->getNumCols();
    137   int numRows = clpSolver->getNumRows();
    138 
    139   // - get the right hand side of the rows
    140   const double * rhs = clpSolver->getRightHandSide();
    141 
    142   // - find the integer variables
    143   bool * varClassInt = new bool[numCols];
    144   int numInt = 0;
    145   for(int i=0; i<numCols; i++)
    146     {
    147       if(clpSolver->isContinuous(i))
    148         varClassInt[i] = 0;
    149       else
    150         {
    151           varClassInt[i] = 1;
    152           numInt++;
    153         }
    154     }
    155  
    156   // -Get the rows sense
    157   const char * rowSense;
    158   rowSense = clpSolver->getRowSense();
    159  
    160   // -Get the objective coefficients
    161   const double *objCoefficients = clpSolver->getObjCoefficients();
    162   double *originalObjCoeff = new double [numCols];
    163   for(int i=0; i<numCols;i++)
    164     originalObjCoeff[i] = objCoefficients[i];
    165  
    166   // -Get the matrix of the problem
    167   double ** matrix = new double * [numRows];
    168   for(int i = 0; i < numRows; i++)
    169     {
    170       matrix[i] = new double[numCols];
    171       for(int j = 0; j < numCols; j++)
    172         matrix[i][j] = 0;
    173     }
    174   const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
    175   const double * matrixElements = matrixByRow->getElements();
    176   const int * matrixIndices = matrixByRow->getIndices();
    177   const int * matrixStarts = matrixByRow->getVectorStarts();
    178   for(int j = 0; j < numRows; j++)
    179     {
    180       for(int i = matrixStarts[j]; i < matrixStarts[j+1]; i++)
    181         {
    182           matrix[j][matrixIndices[i]] = matrixElements[i];
    183         }
    184     }
    185  
    186   // The newObj is the randomly perturbed constraint used to find new
    187   // corner points
    188   double * newObj = new double [numCols];
    189  
    190   // Set the random seed
    191   srand ( time(NULL) +1);
    192   int randNum;
    193  
    194   // We're going to add a new row to the LP formulation
    195   // after finding each new solution.
    196   // Adding a new row requires the new elements and the new indices.
    197   // The elements are original objective function coefficients.
    198   // The indicies are the (dense) columns indices stored in addRowIndex.
    199   // The rhs is the value of the new solution stored in solutionValue.
    200   int * addRowIndex = new int[numCols];
    201   for(int i=0;i<numCols;i++)
    202     addRowIndex[i]=i;
    203 
    204   // The number of feasible solutions found by the PF heuristic.
    205   // This controls the return code of the solution() method.   
    206   int numFeasibles =0;
    207  
    208   // Shuffle the rows
    209   int * index = new int [numRows];
    210   for(int i=0; i<numRows; i++)
    211     index[i]=i;
    212   for(int i=0; i<numRows; i++)
    213     {
    214       int temp = index[i];
    215       int randNumTemp = i+(rand()%(numRows-i));
    216       index[i]=index[randNumTemp];
    217       index[randNumTemp]=temp;
    218     }
    219  
    220   // In the clpSolution struct, we store a lot of column solutions.
    221   // For each perturb objective, we store the solution from each
    222   // iteration of the LP solve.
    223   // For each perturb objective, we look at the collection of
    224   // solutions to do something extremly intelligent :-)
    225   // We could (and should..and will :-) wipe out the block of
    226   // solutions when we're done with them. But for now, we just move on
    227   // and store the next block of solutions for the next (perturbed)
    228   // objective.
    229   // The variable startIndex tells us where the new block begins.
    230   int startIndex = 0;
    231 
    232   // At most "fixThreshold" number of integer variables can be unsatisfied
    233   // for calling smallBranchAndBound().
    234   // The PF Heuristic only fixes fixThreshold number of variables to
    235   // their integer values. Not more. Not less. The reason is to give
    236   // the smallBB some opportunity to find better solutions. If we fix
    237   // everything it might be too many (leading the heuristic to come up
    238   // with infeasibility rather than a useful result).
    239   // (This is an important paramater. And it is dynamically set.)
    240   double fixThreshold;
    241   /*     
    242           if(numInt > 400)
    243           fixThreshold = 17*sqrt(numInt);
    244           if(numInt<=400 && numInt>100)
    245           fixThreshold = 5*sqrt(numInt);
    246           if(numInt<=100)
    247           fixThreshold = 4*sqrt(numInt);
    248   */
    249   // Initialize fixThreshold based on the number of integer
    250   // variables
    251   if(numInt<=100)
    252     fixThreshold = .35*numInt;
    253   if(numInt>100 && numInt<1000)
    254     fixThreshold = .85*numInt;
    255   if(numInt>=1000)
    256     fixThreshold = .1*numInt;
    257 
    258   // Whenever the dynamic system for changing fixThreshold
    259   // kicks in, it changes the parameter by the
    260   // fixThresholdChange amount. 
    261   // (The 25% should be member data and tuned. Another paper!)
    262   double fixThresholdChange = 0.25 * fixThreshold;
    263 
    264   // maxNode is the maximum number of nodes we allow smallBB to
    265   // search. It's initialized to 400 and changed dynamically.
    266   // The 400 should be member data, if we become virtuous.
    267   int maxNode = 400;
    268 
    269   // We control the decision to change maxNode through the boolean
    270   // variable  changeMaxNode. The boolean variable is initialized to
    271   // true and gets set to false under a condition (and is never true
    272   // again.)
    273   // It's flipped off and stays off (in the current incarnation of PF)
    274   bool changeMaxNode = 1;
    275 
    276   // The sumReturnCode is used for the dynamic system that sets
    277   // fixThreshold and changeMaxNode.
    278   //
    279   // We track what's happening in sumReturnCode. There are 8 switches.
    280   // The first 5 switches corresponds to a return code for smallBB.
    281   //
    282   // We want to know how many times we consecutively get the same
    283   // return code.
    284   //
    285   // If "good" return codes are happening often enough, we're happy.
    286   //
    287   // If a "bad" returncodes happen consecutively, we want to
    288   // change something.
    289   //
    290   // The switch 5 is the number of times PF didn't call smallBB
    291   // becuase the number of integer variables that took integer values
    292   // was less than fixThreshold.
    293   //
    294   // The swicth 6 was added for a brilliant idea...to be announced
    295   // later (another paper!)
    296   //
    297   // The switch 7 is the one that changes the max node. Read the
    298   // code. (Todo: Verbalize the brilliant idea for the masses.)
    299   //
    300   int sumReturnCode[8];
    301   /*
    302     sumReturnCode[0] ~ -1 --> problem too big for smallBB
    303     sumReturnCode[1] ~ 0  --> smallBB not finshed and no soln
    304     sumReturnCode[2] ~ 1  --> smallBB not finshed and there is a soln
    305     sumReturnCode[3] ~ 2  --> smallBB finished and no soln
    306     sumReturnCode[4] ~ 3  --> smallBB finished and there is a soln
    307     sumReturnCode[5] ~ didn't call smallBranchAndBound too few to fix
    308     sumReturnCode[6] ~ didn't call smallBranchAndBound too many unsatisfied
    309     sumReturnCode[7] ~ the same as sumReturnCode[1] but becomes zero just if the returnCode is not 0
    310   */
    311  
    312   for(int i=0; i<8; i++)
    313     sumReturnCode[i]=0;
    314   int * colIndex = new int[numCols];
    315   for(int i=0; i<numCols; i++)
    316     colIndex[i]=i;
    317   double cutoff = COIN_DBL_MAX;
    318   bool didMiniBB;
    319 
    320   // Main loop
    321   for(int i=0; i<numRows; i++)
    322     {
    323       // track the number of mini-bb for the dynamic threshold setting
    324       didMiniBB=0;
    325 
    326       for(int k=startIndex; k<solutions.numberSolutions; k++)
    327         //if the point has 0 unsatisfied variables; make sure it is
    328         //feasible. Check integer feasiblity and constraints.
    329         if(solutions.numberUnsatisfied[k] == 0)
    330           {
    331             double feasibility=1;
    332             //check integer feasibility
    333             for(int icol=0; icol<numCols; icol++)
    334               {
    335                 double closest = floor(solutions.solution[k][icol]+0.5);
    336                 if(varClassInt[icol]&&(fabs(solutions.solution[k][icol]-closest)>1e-6))
    337                   {
    338                     feasibility = 0;
    339                     break;     
    340                     }
    341               }
    342             //check if the solution satisfies the constraints
    343             for(int irow = 0; irow < numRows; irow++)
    344               {
    345                 double lhs = 0;
    346                 for(int j = 0; j <numCols; j++)
    347                   lhs += matrix[irow][j] * solutions.solution[k][j];
    348                 if(rowSense[irow] == 'L' && lhs > rhs[irow] + 1e-6)
    349                   {
    350                     feasibility = 0;
    351                     break;
    352                   }
    353                 if(rowSense[irow] == 'G' && lhs < rhs[irow] - 1e-6)
    354                   {
    355                     feasibility = 0;
    356                     break;
    357                   }
    358                 if(rowSense[irow] == 'E' && (lhs - rhs[irow] > 1e-6 || lhs - rhs[irow] < -1e-6))
    359                   {
    360                     feasibility = 0;
    361                     break;
    362                   }
    363               }
    364            
    365             //if feasible, find the objective value and set the cutoff
    366             // for the smallBB and add a new constraint to the LP
    367             // (and update the best solution found so far for the
    368             // return arguments)
    369             if(feasibility)
    370               {
    371                 double objectiveValue=0;
    372                 for(int j=0; j<numCols; j++)
    373                   objectiveValue+=solutions.solution[k][j]*originalObjCoeff[j];
    374                 cutoff=objectiveValue;
    375                 clpSolver->addRow(numCols, addRowIndex, originalObjCoeff,-COIN_DBL_MAX, cutoff);
    376 
    377                 // Todo: pick up the best solution in the block (not
    378                 // the last).
    379                 solutionValue = objectiveValue;
    380                 for(int m=0; m<numCols; m++)
    381                   betterSolution[m] = solutions.solution[k][m];
    382                 numFeasibles++;
    383               }
    384           }
    385      
    386       // Go through the block of solution and decide if to call smallBB
    387       for(int k=startIndex; k<solutions.numberSolutions; k++)
    388         {
    389           if(solutions.numberUnsatisfied[k] <= fixThreshold)
    390             {
    391               // get new copy
    392               OsiSolverInterface * newSolver;
    393               newSolver = new OsiClpSolverInterface(*clpSolver);
    394               newSolver->setObjSense(1);
    395               newSolver->setObjective(originalObjCoeff);
    396               int numberColumns = newSolver->getNumCols();
    397               int numFixed=0;
    398 
    399               // Fix the first fixThreshold number of integer vars
    400               // that are satisfied
    401               for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++)
    402                 {
    403                   if (newSolver->isInteger(iColumn))
    404                     {
    405                       double value = solutions.solution[k][iColumn];
    406                       double intValue = floor(value+0.5);
    407                       if (fabs(value-intValue)<1.0e-5)
    408                         {
    409                           newSolver->setColLower(iColumn,intValue);
    410                           newSolver->setColUpper(iColumn,intValue);
    411                           numFixed++;
    412                           if(numFixed>numInt-fixThreshold)
    413                             break;
    414                         }
    415                     }                 
    416                 }
    417               printf("numFixed: %d\n", numFixed);
    418               printf("fixThreshold: %f\n", fixThreshold);
    419               printf("numInt: %d\n", numInt);
    420               double *newSolution = new double[numCols];
    421               double newSolutionValue;
    422              
    423               // Call smallBB on the modified problem
    424               int returnCode = smallBranchAndBound(newSolver, maxNode, newSolution,
    425                                                    newSolutionValue, cutoff, "mini");
    426 
    427               // If smallBB found a solution, update the better
    428               // solution and solutionValue (we gave smallBB our
    429               // cutoff, so it only finds improving solutions)
    430               if(returnCode == 1 || returnCode == 3)
    431                 {
    432                   numFeasibles ++;
    433                   solutionValue = newSolutionValue;
    434                   for(int m=0; m<numCols; m++)
    435                     betterSolution[m] = newSolution[m];
    436                   printf("cutoff: %f\n", newSolutionValue);
    437                   printf("time: %.2lf\n", CoinCpuTime()-start);
    438                 }
    439               didMiniBB=1;
    440               printf("returnCode: %d\n", returnCode);
    441              
    442               //Update sumReturnCode array
    443               for(int iRC=0; iRC<6; iRC++)
    444                 {
    445                   if(iRC == returnCode+1)
    446                     sumReturnCode[iRC]++;
    447                   else
    448                     sumReturnCode[iRC]=0;                         
    449                 }
    450               if(returnCode!=0)
    451                 sumReturnCode[7]=0;
    452               else
    453                 sumReturnCode[7]++;
    454               if(returnCode==1 || returnCode==3)
    455                 {
    456                   cutoff = newSolutionValue;
    457                   clpSolver->addRow(numCols, addRowIndex, originalObjCoeff,-COIN_DBL_MAX, cutoff);
    458                   printf("******************\n\n*****************\n");                 
    459                 }
    460               break;
    461             }
    462         }
    463      
    464       if(!didMiniBB && solutions.numberSolutions - startIndex > 0)
    465         {
    466           sumReturnCode[5]++;
    467           for(int iRC=0; iRC<5; iRC++)
    468             sumReturnCode[iRC]=0;
    469         }
    470 
    471       //Change "fixThreshold" if needed
    472       // using the data we've recorded in sumReturnCode
    473       if(sumReturnCode[1]>=3)
    474         fixThreshold -= fixThresholdChange;
    475       if(sumReturnCode[7]>=3 && changeMaxNode)
    476         {
    477           maxNode *= 5;
    478           changeMaxNode=0;
    479         }
    480       if(sumReturnCode[3]>=3 && fixThreshold < 0.95 * numInt)
    481         fixThreshold += fixThresholdChange;
    482       if(sumReturnCode[5]>=4)
    483             fixThreshold += fixThresholdChange;
    484       if(sumReturnCode[0]>3)
    485         fixThreshold -= fixThresholdChange;
    486      
    487       startIndex = solutions.numberSolutions;
    488      
    489       //Check if the maximum iterations limit is reached
    490       // rlh: Ask John how this is working with the change to trustedUserPtr.
    491       if(solutions.numberSolutions > 20000)
    492         break;
    493 
    494       // The first time in this loop PF solves orig LP.
    495      
    496       //Generate the random objective function
    497       randNum = rand()%10 + 1;
    498       randNum = fmod(randNum, 2);
    499       for(int j=0; j<numCols; j++)
    500         {
    501           if(randNum == 1)             
    502             if(fabs(matrix[index[i]][j]) < 1e-6)
    503               newObj[j] = 0.1;
    504             else
    505               newObj[j] = matrix[index[i]][j] * 1.1;
    506           else
    507             if(fabs(matrix[index[i]][j]) < 1e-6)
    508               newObj[j] = -0.1;
    509             else
    510               newObj[j] = matrix[index[i]][j] * 0.9;
    511         }
    512       clpSolver->setObjective(newObj);
    513       if(rowSense[i] == 'L')
    514         clpSolver->setObjSense(-1);
    515       else
    516         // Todo #1: We don't need to solve the LPs to optimality.
    517         // We just need corner points.
    518         // There's a problem in stopping Clp that needs to be looked
    519         // into. So for now, we solve optimality.
    520         clpSolver->setObjSense(1);
    521       //          simplex->setMaximumIterations(100);
    522       clpSolver->getModelPtr()->primal(1);
    523       //          simplex->setMaximumIterations(100000);
    524       printf("cutoff: %f\n", cutoff);
    525       printf("time: %.2f\n", CoinCpuTime()-start);
    526       for(int iRC=0; iRC<8; iRC++)
    527           printf("%d ",sumReturnCode[iRC]);
    528       printf("\nfixThreshold: %f\n", fixThreshold);
    529       printf("numInt: %d\n", numInt);
    530       printf("\n---------------------------------------------------------------- %d\n", i);
    531 
    532       //temp:
    533       if(i>3) break;
    534 
    535     }
    536 
    537   printf("Best Feasible Found: %f\n", cutoff);
    538   printf("Total time: %.2f\n", CoinCpuTime()-start);
    539  
    540   if (numFeasibles == 0) {
    541     return 0;
    542   }
    543  
    544 
    545 
    546   // We found something better
    547   std::cout << "See you soon! You're leaving the Pivot-and-Fix Heuristic" << std::endl;
    548  std::cout << std::endl;
    549 
    550   return 1;
     110    assert (clpSolverOriginal);
     111    OsiClpSolverInterface *clpSolver(clpSolverOriginal);
     112
     113    ClpSimplex * simplex = clpSolver->getModelPtr();
     114
     115    // Initialize the structure holding the solutions
     116    clpSolution solutions;
     117    // Set typeStruct field of ClpTrustedData struct to one.
     118    // This tells Clp it's "Mahdi!"
     119    ClpTrustedData trustedSolutions;
     120    trustedSolutions.typeStruct = 1;
     121    trustedSolutions.data = &solutions;
     122    solutions.numberSolutions = 0;
     123    solutions.maximumSolutions = 0;
     124    solutions.numberColumns = simplex->numberColumns();
     125    solutions.solution = NULL;
     126    solutions.numberUnsatisfied = NULL;
     127    simplex->setTrustedUserPointer(&trustedSolutions);
     128
     129    // Solve from all slack to get some points
     130    simplex->allSlackBasis();
     131    simplex->primal();
     132
     133    // -------------------------------------------------
     134    // Get the problem information
     135    // - get the number of cols and rows
     136    int numCols = clpSolver->getNumCols();
     137    int numRows = clpSolver->getNumRows();
     138
     139    // - get the right hand side of the rows
     140    const double * rhs = clpSolver->getRightHandSide();
     141
     142    // - find the integer variables
     143    bool * varClassInt = new bool[numCols];
     144    int numInt = 0;
     145    for (int i = 0; i < numCols; i++) {
     146        if (clpSolver->isContinuous(i))
     147            varClassInt[i] = 0;
     148        else {
     149            varClassInt[i] = 1;
     150            numInt++;
     151        }
     152    }
     153
     154    // -Get the rows sense
     155    const char * rowSense;
     156    rowSense = clpSolver->getRowSense();
     157
     158    // -Get the objective coefficients
     159    const double *objCoefficients = clpSolver->getObjCoefficients();
     160    double *originalObjCoeff = new double [numCols];
     161    for (int i = 0; i < numCols; i++)
     162        originalObjCoeff[i] = objCoefficients[i];
     163
     164    // -Get the matrix of the problem
     165    double ** matrix = new double * [numRows];
     166    for (int i = 0; i < numRows; i++) {
     167        matrix[i] = new double[numCols];
     168        for (int j = 0; j < numCols; j++)
     169            matrix[i][j] = 0;
     170    }
     171    const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
     172    const double * matrixElements = matrixByRow->getElements();
     173    const int * matrixIndices = matrixByRow->getIndices();
     174    const int * matrixStarts = matrixByRow->getVectorStarts();
     175    for (int j = 0; j < numRows; j++) {
     176        for (int i = matrixStarts[j]; i < matrixStarts[j+1]; i++) {
     177            matrix[j][matrixIndices[i]] = matrixElements[i];
     178        }
     179    }
     180
     181    // The newObj is the randomly perturbed constraint used to find new
     182    // corner points
     183    double * newObj = new double [numCols];
     184
     185    // Set the random seed
     186    srand ( time(NULL) + 1);
     187    int randNum;
     188
     189    // We're going to add a new row to the LP formulation
     190    // after finding each new solution.
     191    // Adding a new row requires the new elements and the new indices.
     192    // The elements are original objective function coefficients.
     193    // The indicies are the (dense) columns indices stored in addRowIndex.
     194    // The rhs is the value of the new solution stored in solutionValue.
     195    int * addRowIndex = new int[numCols];
     196    for (int i = 0; i < numCols; i++)
     197        addRowIndex[i] = i;
     198
     199    // The number of feasible solutions found by the PF heuristic.
     200    // This controls the return code of the solution() method.
     201    int numFeasibles = 0;
     202
     203    // Shuffle the rows
     204    int * index = new int [numRows];
     205    for (int i = 0; i < numRows; i++)
     206        index[i] = i;
     207    for (int i = 0; i < numRows; i++) {
     208        int temp = index[i];
     209        int randNumTemp = i + (rand() % (numRows - i));
     210        index[i] = index[randNumTemp];
     211        index[randNumTemp] = temp;
     212    }
     213
     214    // In the clpSolution struct, we store a lot of column solutions.
     215    // For each perturb objective, we store the solution from each
     216    // iteration of the LP solve.
     217    // For each perturb objective, we look at the collection of
     218    // solutions to do something extremly intelligent :-)
     219    // We could (and should..and will :-) wipe out the block of
     220    // solutions when we're done with them. But for now, we just move on
     221    // and store the next block of solutions for the next (perturbed)
     222    // objective.
     223    // The variable startIndex tells us where the new block begins.
     224    int startIndex = 0;
     225
     226    // At most "fixThreshold" number of integer variables can be unsatisfied
     227    // for calling smallBranchAndBound().
     228    // The PF Heuristic only fixes fixThreshold number of variables to
     229    // their integer values. Not more. Not less. The reason is to give
     230    // the smallBB some opportunity to find better solutions. If we fix
     231    // everything it might be too many (leading the heuristic to come up
     232    // with infeasibility rather than a useful result).
     233    // (This is an important paramater. And it is dynamically set.)
     234    double fixThreshold;
     235    /*
     236        if(numInt > 400)
     237        fixThreshold = 17*sqrt(numInt);
     238        if(numInt<=400 && numInt>100)
     239        fixThreshold = 5*sqrt(numInt);
     240        if(numInt<=100)
     241        fixThreshold = 4*sqrt(numInt);
     242    */
     243    // Initialize fixThreshold based on the number of integer
     244    // variables
     245    if (numInt <= 100)
     246        fixThreshold = .35 * numInt;
     247    if (numInt > 100 && numInt < 1000)
     248        fixThreshold = .85 * numInt;
     249    if (numInt >= 1000)
     250        fixThreshold = .1 * numInt;
     251
     252    // Whenever the dynamic system for changing fixThreshold
     253    // kicks in, it changes the parameter by the
     254    // fixThresholdChange amount.
     255    // (The 25% should be member data and tuned. Another paper!)
     256    double fixThresholdChange = 0.25 * fixThreshold;
     257
     258    // maxNode is the maximum number of nodes we allow smallBB to
     259    // search. It's initialized to 400 and changed dynamically.
     260    // The 400 should be member data, if we become virtuous.
     261    int maxNode = 400;
     262
     263    // We control the decision to change maxNode through the boolean
     264    // variable  changeMaxNode. The boolean variable is initialized to
     265    // true and gets set to false under a condition (and is never true
     266    // again.)
     267    // It's flipped off and stays off (in the current incarnation of PF)
     268    bool changeMaxNode = 1;
     269
     270    // The sumReturnCode is used for the dynamic system that sets
     271    // fixThreshold and changeMaxNode.
     272    //
     273    // We track what's happening in sumReturnCode. There are 8 switches.
     274    // The first 5 switches corresponds to a return code for smallBB.
     275    //
     276    // We want to know how many times we consecutively get the same
     277    // return code.
     278    //
     279    // If "good" return codes are happening often enough, we're happy.
     280    //
     281    // If a "bad" returncodes happen consecutively, we want to
     282    // change something.
     283    //
     284    // The switch 5 is the number of times PF didn't call smallBB
     285    // becuase the number of integer variables that took integer values
     286    // was less than fixThreshold.
     287    //
     288    // The swicth 6 was added for a brilliant idea...to be announced
     289    // later (another paper!)
     290    //
     291    // The switch 7 is the one that changes the max node. Read the
     292    // code. (Todo: Verbalize the brilliant idea for the masses.)
     293    //
     294    int sumReturnCode[8];
     295    /*
     296      sumReturnCode[0] ~ -1 --> problem too big for smallBB
     297      sumReturnCode[1] ~ 0  --> smallBB not finshed and no soln
     298      sumReturnCode[2] ~ 1  --> smallBB not finshed and there is a soln
     299      sumReturnCode[3] ~ 2  --> smallBB finished and no soln
     300      sumReturnCode[4] ~ 3  --> smallBB finished and there is a soln
     301      sumReturnCode[5] ~ didn't call smallBranchAndBound too few to fix
     302      sumReturnCode[6] ~ didn't call smallBranchAndBound too many unsatisfied
     303      sumReturnCode[7] ~ the same as sumReturnCode[1] but becomes zero just if the returnCode is not 0
     304    */
     305
     306    for (int i = 0; i < 8; i++)
     307        sumReturnCode[i] = 0;
     308    int * colIndex = new int[numCols];
     309    for (int i = 0; i < numCols; i++)
     310        colIndex[i] = i;
     311    double cutoff = COIN_DBL_MAX;
     312    bool didMiniBB;
     313
     314    // Main loop
     315    for (int i = 0; i < numRows; i++) {
     316        // track the number of mini-bb for the dynamic threshold setting
     317        didMiniBB = 0;
     318
     319        for (int k = startIndex; k < solutions.numberSolutions; k++)
     320            //if the point has 0 unsatisfied variables; make sure it is
     321            //feasible. Check integer feasiblity and constraints.
     322            if (solutions.numberUnsatisfied[k] == 0) {
     323                double feasibility = 1;
     324                //check integer feasibility
     325                for (int icol = 0; icol < numCols; icol++) {
     326                    double closest = floor(solutions.solution[k][icol] + 0.5);
     327                    if (varClassInt[icol] && (fabs(solutions.solution[k][icol] - closest) > 1e-6)) {
     328                        feasibility = 0;
     329                        break;
     330                    }
     331                }
     332                //check if the solution satisfies the constraints
     333                for (int irow = 0; irow < numRows; irow++) {
     334                    double lhs = 0;
     335                    for (int j = 0; j < numCols; j++)
     336                        lhs += matrix[irow][j] * solutions.solution[k][j];
     337                    if (rowSense[irow] == 'L' && lhs > rhs[irow] + 1e-6) {
     338                        feasibility = 0;
     339                        break;
     340                    }
     341                    if (rowSense[irow] == 'G' && lhs < rhs[irow] - 1e-6) {
     342                        feasibility = 0;
     343                        break;
     344                    }
     345                    if (rowSense[irow] == 'E' && (lhs - rhs[irow] > 1e-6 || lhs - rhs[irow] < -1e-6)) {
     346                        feasibility = 0;
     347                        break;
     348                    }
     349                }
     350
     351                //if feasible, find the objective value and set the cutoff
     352                // for the smallBB and add a new constraint to the LP
     353                // (and update the best solution found so far for the
     354                // return arguments)
     355                if (feasibility) {
     356                    double objectiveValue = 0;
     357                    for (int j = 0; j < numCols; j++)
     358                        objectiveValue += solutions.solution[k][j] * originalObjCoeff[j];
     359                    cutoff = objectiveValue;
     360                    clpSolver->addRow(numCols, addRowIndex, originalObjCoeff, -COIN_DBL_MAX, cutoff);
     361
     362                    // Todo: pick up the best solution in the block (not
     363                    // the last).
     364                    solutionValue = objectiveValue;
     365                    for (int m = 0; m < numCols; m++)
     366                        betterSolution[m] = solutions.solution[k][m];
     367                    numFeasibles++;
     368                }
     369            }
     370
     371        // Go through the block of solution and decide if to call smallBB
     372        for (int k = startIndex; k < solutions.numberSolutions; k++) {
     373            if (solutions.numberUnsatisfied[k] <= fixThreshold) {
     374                // get new copy
     375                OsiSolverInterface * newSolver;
     376                newSolver = new OsiClpSolverInterface(*clpSolver);
     377                newSolver->setObjSense(1);
     378                newSolver->setObjective(originalObjCoeff);
     379                int numberColumns = newSolver->getNumCols();
     380                int numFixed = 0;
     381
     382                // Fix the first fixThreshold number of integer vars
     383                // that are satisfied
     384                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
     385                    if (newSolver->isInteger(iColumn)) {
     386                        double value = solutions.solution[k][iColumn];
     387                        double intValue = floor(value + 0.5);
     388                        if (fabs(value - intValue) < 1.0e-5) {
     389                            newSolver->setColLower(iColumn, intValue);
     390                            newSolver->setColUpper(iColumn, intValue);
     391                            numFixed++;
     392                            if (numFixed > numInt - fixThreshold)
     393                                break;
     394                        }
     395                    }
     396                }
     397                printf("numFixed: %d\n", numFixed);
     398                printf("fixThreshold: %f\n", fixThreshold);
     399                printf("numInt: %d\n", numInt);
     400                double *newSolution = new double[numCols];
     401                double newSolutionValue;
     402
     403                // Call smallBB on the modified problem
     404                int returnCode = smallBranchAndBound(newSolver, maxNode, newSolution,
     405                                                     newSolutionValue, cutoff, "mini");
     406
     407                // If smallBB found a solution, update the better
     408                // solution and solutionValue (we gave smallBB our
     409                // cutoff, so it only finds improving solutions)
     410                if (returnCode == 1 || returnCode == 3) {
     411                    numFeasibles ++;
     412                    solutionValue = newSolutionValue;
     413                    for (int m = 0; m < numCols; m++)
     414                        betterSolution[m] = newSolution[m];
     415                    printf("cutoff: %f\n", newSolutionValue);
     416                    printf("time: %.2lf\n", CoinCpuTime() - start);
     417                }
     418                didMiniBB = 1;
     419                printf("returnCode: %d\n", returnCode);
     420
     421                //Update sumReturnCode array
     422                for (int iRC = 0; iRC < 6; iRC++) {
     423                    if (iRC == returnCode + 1)
     424                        sumReturnCode[iRC]++;
     425                    else
     426                        sumReturnCode[iRC] = 0;
     427                }
     428                if (returnCode != 0)
     429                    sumReturnCode[7] = 0;
     430                else
     431                    sumReturnCode[7]++;
     432                if (returnCode == 1 || returnCode == 3) {
     433                    cutoff = newSolutionValue;
     434                    clpSolver->addRow(numCols, addRowIndex, originalObjCoeff, -COIN_DBL_MAX, cutoff);
     435                    printf("******************\n\n*****************\n");
     436                }
     437                break;
     438            }
     439        }
     440
     441        if (!didMiniBB && solutions.numberSolutions - startIndex > 0) {
     442            sumReturnCode[5]++;
     443            for (int iRC = 0; iRC < 5; iRC++)
     444                sumReturnCode[iRC] = 0;
     445        }
     446
     447        //Change "fixThreshold" if needed
     448        // using the data we've recorded in sumReturnCode
     449        if (sumReturnCode[1] >= 3)
     450            fixThreshold -= fixThresholdChange;
     451        if (sumReturnCode[7] >= 3 && changeMaxNode) {
     452            maxNode *= 5;
     453            changeMaxNode = 0;
     454        }
     455        if (sumReturnCode[3] >= 3 && fixThreshold < 0.95 * numInt)
     456            fixThreshold += fixThresholdChange;
     457        if (sumReturnCode[5] >= 4)
     458            fixThreshold += fixThresholdChange;
     459        if (sumReturnCode[0] > 3)
     460            fixThreshold -= fixThresholdChange;
     461
     462        startIndex = solutions.numberSolutions;
     463
     464        //Check if the maximum iterations limit is reached
     465        // rlh: Ask John how this is working with the change to trustedUserPtr.
     466        if (solutions.numberSolutions > 20000)
     467            break;
     468
     469        // The first time in this loop PF solves orig LP.
     470
     471        //Generate the random objective function
     472        randNum = rand() % 10 + 1;
     473        randNum = fmod(randNum, 2);
     474        for (int j = 0; j < numCols; j++) {
     475            if (randNum == 1)
     476                if (fabs(matrix[index[i]][j]) < 1e-6)
     477                    newObj[j] = 0.1;
     478                else
     479                    newObj[j] = matrix[index[i]][j] * 1.1;
     480            else if (fabs(matrix[index[i]][j]) < 1e-6)
     481                newObj[j] = -0.1;
     482            else
     483                newObj[j] = matrix[index[i]][j] * 0.9;
     484        }
     485        clpSolver->setObjective(newObj);
     486        if (rowSense[i] == 'L')
     487            clpSolver->setObjSense(-1);
     488        else
     489            // Todo #1: We don't need to solve the LPs to optimality.
     490            // We just need corner points.
     491            // There's a problem in stopping Clp that needs to be looked
     492            // into. So for now, we solve optimality.
     493            clpSolver->setObjSense(1);
     494        //        simplex->setMaximumIterations(100);
     495        clpSolver->getModelPtr()->primal(1);
     496        //        simplex->setMaximumIterations(100000);
     497        printf("cutoff: %f\n", cutoff);
     498        printf("time: %.2f\n", CoinCpuTime() - start);
     499        for (int iRC = 0; iRC < 8; iRC++)
     500            printf("%d ", sumReturnCode[iRC]);
     501        printf("\nfixThreshold: %f\n", fixThreshold);
     502        printf("numInt: %d\n", numInt);
     503        printf("\n---------------------------------------------------------------- %d\n", i);
     504
     505        //temp:
     506        if (i > 3) break;
     507
     508    }
     509
     510    printf("Best Feasible Found: %f\n", cutoff);
     511    printf("Total time: %.2f\n", CoinCpuTime() - start);
     512
     513    if (numFeasibles == 0) {
     514        return 0;
     515    }
     516
     517
     518
     519    // We found something better
     520    std::cout << "See you soon! You're leaving the Pivot-and-Fix Heuristic" << std::endl;
     521    std::cout << std::endl;
     522
     523    return 1;
    551524#endif
    552525
    553   return 0;
     526    return 0;
    554527
    555528}
     
    561534void CbcHeuristicPivotAndFix::setModel(CbcModel * )
    562535{
    563   // probably same as resetModel
    564 }
    565 
    566  
     536    // probably same as resetModel
     537}
     538
     539
Note: See TracChangeset for help on using the changeset viewer.