Changeset 1065


Ignore:
Timestamp:
Sep 12, 2008 5:08:57 PM (10 years ago)
Author:
rlh
Message:

adding guts of CbcHeuristicPivotFix?

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcHeuristicPivotAndFix.cpp

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