Changeset 1064


Ignore:
Timestamp:
Sep 11, 2008 4:56:14 PM (10 years ago)
Author:
rlh
Message:

adding the guts of the randomized rounding heuristic

File:
1 edited

Legend:

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

    r1057 r1064  
    99#include <cmath>
    1010#include <cfloat>
     11#include <vector.h>
    1112
    1213#include "OsiSolverInterface.hpp"
     
    1415#include "CbcMessage.hpp"
    1516#include "CbcHeuristicRandRound.hpp"
     17#include "OsiClpSolverInterface.hpp"
     18#include  "CoinTime.hpp"
    1619
    1720// Default Constructor
     
    7477}
    7578/*
    76   Comments needed
     79  Randomized Rounding Heuristic
    7780  Returns 1 if solution, 0 if not */
    7881int
     
    8083                         double * betterSolution)
    8184{
    82 
    83   numCouldRun_++;
    84   int returnCode = 0;
    85   return returnCode;
     85  std::cout << "Lucky you! You're in the Randomized Rounding Heuristic" << std::endl;
     86  // The struct should be moved to member data
     87
     88  typedef struct {
     89    int numberSolutions;
     90    int maximumSolutions;
     91    int numberColumns;
     92    double ** solution;
     93    int * numberUnsatisfied;
     94  } clpSolution;
     95
     96  double start = CoinCpuTime();
     97  numCouldRun_++; // Todo: Ask JJHF what this for.
     98
     99  OsiClpSolverInterface * clpSolver
     100    = dynamic_cast<OsiClpSolverInterface *> (model_->solver());
     101  assert (clpSolver);
     102  ClpSimplex * simplex = clpSolver->getModelPtr();
     103
     104  // Initialize the structure holding the solutions
     105  clpSolution solutions;
     106  // Set typeStruct field of ClpTrustedData struct to one.
     107  // This tells Clp it's "Mahdi!"
     108  ClpTrustedData trustedSolutions;
     109  trustedSolutions.typeStruct = 1;
     110  trustedSolutions.data = &solutions;
     111  solutions.numberSolutions=0;
     112  solutions.maximumSolutions=0;
     113  solutions.numberColumns=simplex->numberColumns();
     114  solutions.solution=NULL;
     115  solutions.numberUnsatisfied=NULL;
     116  simplex->setTrustedUserPointer(&trustedSolutions); // rlh: old was
     117                                              // "userPointer"
     118   // Solve from all slack to get some points
     119  simplex->allSlackBasis();
     120  simplex->primal();
     121
     122  // rlh: Mahdi's code
     123  // -------------------------------------------------
     124  // Get the problem information
     125 
     126  // - Get the number of cols and rows
     127  int numCols = clpSolver->getNumCols();
     128  int numRows = clpSolver->getNumRows();
     129
     130  // - Get the right hand side of the rows
     131  const double * rhs = clpSolver->getRightHandSide();
     132
     133  // Find the integer variables
     134  // rlh: consider using columnType
     135  // One if not continuous (i.e., bin or gen int)
     136  bool * varClassInt = new bool[numCols];
     137  for(int i=0; i<numCols; i++)
     138    {
     139      if(clpSolver->isContinuous(i))
     140        varClassInt[i] = 0;
     141      else
     142        varClassInt[i] = 1;
     143    }
     144
     145  // -Get the rows sense
     146  const char * rowSense;
     147  rowSense = clpSolver->getRowSense();
     148
     149  // -Get the objective coefficients
     150  const double *objCoefficients = clpSolver->getObjCoefficients();
     151  double *originalObjCoeff = new double [numCols];
     152  for(int i=0; i<numCols;i++)
     153    originalObjCoeff[i] = objCoefficients[i];
     154
     155  // -Get the matrix of the problem
     156  double ** matrix = new double * [numRows];
     157  for(int i = 0; i < numRows; i++)
     158    {
     159      matrix[i] = new double[numCols];
     160      for(int j = 0; j < numCols; j++)
     161        matrix[i][j] = 0;
     162    }
     163 
     164  const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
     165  const double * matrixElements = matrixByRow->getElements();
     166  const int * matrixIndices = matrixByRow->getIndices();
     167  const int * matrixStarts = matrixByRow->getVectorStarts();
     168  for(int j = 0; j < numRows; j++)
     169    {
     170      for(int i = matrixStarts[j]; i < matrixStarts[j+1]; i++)
     171        {
     172          matrix[j][matrixIndices[i]] = matrixElements[i];
     173        }
     174    }
     175 
     176  double * newObj = new double [numCols];
     177  srand ( time(NULL) +1);
     178  int randNum;
     179
     180  // Shuffle the rows:
     181  // Put the rows in a random order
     182  // so that the optimal solution is a different corner point than the
     183  // starting point.
     184  int * index = new int [numRows];
     185  for(int i=0; i<numRows; i++)
     186    index[i]=i;
     187  for(int i=0; i<numRows; i++)
     188    {
     189      int temp = index[i];
     190      int randNumTemp = i+(rand()%(numRows-i));
     191      index[i]=index[randNumTemp];
     192      index[randNumTemp]=temp;
     193    }
     194 
     195  // Start finding corner points by iteratively doing the following:
     196  // - find randomly tilted objective
     197  // - solve
     198  for(int i=0; i<numRows; i++)
     199    {
     200      // TODO: that 10,000 could be a param in the member data
     201      if(solutions.numberSolutions  > 10000)
     202        break;
     203      randNum = rand()%10 + 1;
     204      randNum = fmod(randNum, 2);
     205      for(int j=0; j<numCols; j++)
     206        {
     207          // for row i and column j vary the coefficient "a bit"
     208          if(randNum == 1)             
     209            // if the element is zero, then round the coefficient down to 0.1
     210            if(fabs(matrix[index[i]][j]) < 1e-6)
     211              newObj[j] = 0.1;
     212            else
     213              // if the element is nonzero, then increase it "a bit"
     214              newObj[j] = matrix[index[i]][j] * 1.1;
     215          else
     216            // if randnum is 2, then
     217            // if the element is zero, then round the coefficient down
     218            // to NEGATIVE 0.1
     219            if(fabs(matrix[index[i]][j]) < 1e-6)
     220              newObj[j] = -0.1;
     221            else
     222              // if the element is nonzero, then DEcrease it "a bit"
     223              newObj[j] = matrix[index[i]][j] * 0.9;
     224        }
     225      // Use the new "tilted" objective
     226      clpSolver->setObjective(newObj);
     227
     228      // Based on the row sense, we decide whether to max or min
     229      if(rowSense[i] == 'L')
     230        clpSolver->setObjSense(-1);
     231      else
     232        clpSolver->setObjSense(1);
     233
     234      // Solve with primal simplex
     235      clpSolver->getModelPtr()->primal(1);
     236      printf("---------------------------------------------------------------- %d\n", i);
     237    }
     238  // Iteratively do this process until...
     239  // either you reach the max number of corner points (aka 10K)
     240  // or all the rows have been used as an objective.
     241 
     242  // Look at solutions
     243  int numberSolutions = solutions.numberSolutions;
     244  const char * integerInfo = simplex->integerInformation();
     245  const double * columnLower = simplex->columnLower();
     246  const double * columnUpper = simplex->columnUpper();
     247  printf("there are %d solutions\n",numberSolutions);
     248
     249  // Up to here we have all the corner points
     250  // Now we need to do the random walks and roundings
     251 
     252  double ** cornerPoints = new double * [numberSolutions];
     253  for(int j=0;j<numberSolutions;j++)
     254    cornerPoints[j] = solutions.solution[j];
     255 
     256  bool feasibility = 1;
     257  double bestObj = 1e30;
     258  vector< vector <double> > feasibles;
     259  int numFeasibles = 0;
     260 
     261  // Check the feasibility of the corner points
     262  int numCornerPoints = numberSolutions;
     263
     264  rhs = clpSolver->getRightHandSide();
     265  rowSense = clpSolver->getRowSense();
     266 
     267  for(int i=0; i<numCornerPoints; i++)
     268    {
     269      //get the objective value for this this point
     270      double objValue = 0;
     271      for(int k=0; k<numCols; k++)
     272        objValue += cornerPoints[i][k] * originalObjCoeff[k];
     273     
     274      if(objValue < bestObj)
     275        {
     276          feasibility = 1;
     277          for(int j=0; j<numCols; j++)
     278            {
     279              if(varClassInt[j])
     280                {
     281                  double closest=floor(cornerPoints[i][j]+0.5);
     282                  if(fabs(cornerPoints[i][j] - closest)>1e-6)
     283                    {
     284                      feasibility = 0;
     285                      break;
     286                    }
     287                }
     288            }
     289          if(feasibility)
     290              {
     291                for(int irow = 0; irow < numRows; irow++)
     292                  {
     293                    double lhs = 0;
     294                    for(int j = 0; j <numCols; j++)
     295                      {
     296                        lhs += matrix[irow][j] * cornerPoints[i][j];
     297                      }
     298                    if(rowSense[irow] == 'L' && lhs > rhs[irow] + 1e-6)
     299                      {
     300                        feasibility = 0;
     301                        break;
     302                      }
     303                    if(rowSense[irow] == 'G' && lhs < rhs[irow] - 1e-6)
     304                      {
     305                        feasibility = 0;
     306                        break;
     307                      }
     308                    if(rowSense[irow] == 'E' && (lhs - rhs[irow] > 1e-6 || lhs - rhs[irow] < -1e-6))
     309                      {
     310                        feasibility = 0;
     311                        break;
     312                      }
     313                  }
     314              }
     315         
     316          if(feasibility)
     317            {           
     318              numFeasibles++;
     319                feasibles.push_back(vector <double> (numCols));
     320                for(int k=0; k<numCols; k++)
     321                  feasibles[numFeasibles-1][k] = cornerPoints[i][k];
     322                printf("obj: %f\n", objValue);
     323                if(objValue < bestObj)
     324                  bestObj = objValue;
     325            }
     326        }
     327    }   
     328  int numFeasibleCorners;
     329  numFeasibleCorners = numFeasibles;
     330  //find the center of gravity of the corner points as the first random point
     331  double * rp = new double[numCols];
     332  for(int i=0; i<numCols; i++)
     333    {
     334      rp[i] = 0;
     335      for(int j=0; j < numCornerPoints; j++)
     336        {
     337            rp[i] += cornerPoints[j][i];
     338        }
     339      rp[i] = rp[i]/numCornerPoints;
     340    }
     341 
     342  //-------------------------------------------
     343  //main loop:
     344  // -generate the next random point
     345  // -round the random point
     346  // -check the feasibility of the random point
     347  //-------------------------------------------
     348 
     349  srand ( time(NULL) +1);
     350  int numRandomPoints = 0;
     351  while(numRandomPoints < 50000)
     352    {
     353      numRandomPoints++;
     354      //generate the next random point
     355      int randomIndex = rand()%numCornerPoints;
     356      double random = rand();
     357      random = random/RAND_MAX;
     358      for(int i=0; i<numCols; i++)
     359        {
     360          rp[i] = (random * (cornerPoints[randomIndex][i] - rp[i])) + rp[i];
     361        }
     362     
     363      //CRISP ROUNDING
     364      //round the random point just generated
     365      double * roundRp = new double[numCols];
     366      for(int i=0; i<numCols; i++)
     367        {
     368          roundRp[i] = rp[i];
     369          if(varClassInt[i])
     370            {
     371              if(rp[i]>=0)
     372                {                               
     373                  if(fmod(rp[i],1) > 0.5)
     374                    roundRp[i] = floor(rp[i]) + 1;
     375                  else
     376                    roundRp[i] = floor(rp[i]);
     377                }
     378              else
     379                {
     380                  if(abs(fmod(rp[i],1)) > 0.5)
     381                    roundRp[i] = floor(rp[i]);
     382                  else
     383                    roundRp[i] = floor(rp[i])+1;
     384                 
     385                }
     386            }
     387        }
     388     
     389     
     390      //SOFT ROUNDING
     391      // Look at original files for the "how to" on soft rounding
     392     
     393     
     394      //check the feasibility of the rounded random point
     395      // -Check the feasibility
     396      // -Get the rows sense
     397      rowSense = clpSolver->getRowSense();
     398      rhs = clpSolver->getRightHandSide();
     399     
     400      //get the objective value for this feasible point
     401      double objValue = 0;
     402      for(int i=0; i<numCols; i++)
     403        objValue += roundRp[i] * originalObjCoeff[i];
     404     
     405      if(objValue<bestObj)
     406        {
     407          feasibility = 1;
     408          for(int i = 0; i < numRows; i++)
     409            {
     410              double lhs = 0;
     411              for(int j = 0; j <numCols; j++)
     412                {
     413                  lhs += matrix[i][j] * roundRp[j];
     414                }
     415              if(rowSense[i] == 'L' && lhs > rhs[i] + 1e-6)
     416                {
     417                  feasibility = 0;
     418                  break;
     419                }
     420              if(rowSense[i] == 'G' && lhs < rhs[i] - 1e-6)
     421                {
     422                  feasibility = 0;
     423                  break;
     424                }
     425              if(rowSense[i] == 'E' && (lhs - rhs[i] > 1e-6 || lhs - rhs[i] < -1e-6))
     426                {
     427                  feasibility = 0;
     428                  break;
     429                }
     430            }
     431          if(feasibility)
     432            {
     433              printf("Feasible Found!!\n");
     434              printf("%.2f\n", CoinCpuTime()-start);
     435              numFeasibles++;
     436              feasibles.push_back(vector <double> (numCols));
     437              for(int i=0; i<numCols; i++)
     438                feasibles[numFeasibles-1][i] = roundRp[i];
     439              printf("obj: %f\n", objValue);
     440              if(objValue < bestObj)
     441                bestObj = objValue;
     442            }
     443        }
     444    }
     445  printf("Number of Feasible Corners: %d\n", numFeasibleCorners);
     446  printf("Number of Feasibles Found: %d\n", numFeasibles);
     447  if(numFeasibles > 0)
     448    printf("Best Objective: %f\n", bestObj);
     449  printf("time: %.2f\n",CoinCpuTime()-start);
     450 
     451  if (numFeasibles == 0) return 0;
     452 
     453  // We found something better
     454  solutionValue = bestObj;
     455  for (int k = 0; k<numCols; k++){
     456    betterSolution[k] =  feasibles[numFeasibles-1][k];
     457  } 
     458  std::cout << "See you soon! You're leaving the Randomized Rounding Heuristic" << std::endl;
     459  return 1;
     460
    86461}
    87462// update model
Note: See TracChangeset for help on using the changeset viewer.