source: branches/primalheur/Cbc/src/CbcHeuristicRandRound.cpp @ 1109

Last change on this file since 1109 was 1109, checked in by mahdi, 13 years ago

test

File size: 14.8 KB
Line 
1//heeeeeeeeeeeeey
2// Copyright (C) 2008, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12#include <vector>
13
14#include "CoinHelperFunctions.hpp"
15#include "OsiSolverInterface.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcHeuristicRandRound.hpp"
19#include "OsiClpSolverInterface.hpp"
20#include  "CoinTime.hpp"
21
22static inline int intRand(const int range)
23{
24  return (int) floor(CoinDrand48() * range);
25}
26
27// Default Constructor
28CbcHeuristicRandRound::CbcHeuristicRandRound() 
29  :CbcHeuristic()
30{
31}
32
33// Constructor with model - assumed before cuts
34CbcHeuristicRandRound::CbcHeuristicRandRound(CbcModel & model)
35  :CbcHeuristic(model)
36{
37  model_=&model;
38  setWhen(1);
39}
40
41// Destructor
42CbcHeuristicRandRound::~CbcHeuristicRandRound ()
43{
44}
45
46// Clone
47CbcHeuristic *
48CbcHeuristicRandRound::clone() const
49{
50  return new CbcHeuristicRandRound(*this);
51}
52
53// Create C++ lines to get to current state
54void 
55CbcHeuristicRandRound::generateCpp( FILE * fp) 
56{
57  CbcHeuristicRandRound other;
58  fprintf(fp,"0#include \"CbcHeuristicRandRound.hpp\"\n");
59  fprintf(fp,"3  CbcHeuristicRandRound heuristicPFX(*cbcModel);\n");
60  CbcHeuristic::generateCpp(fp,"heuristicPFX");
61  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicPFX);\n");
62}
63
64// Copy constructor
65CbcHeuristicRandRound::CbcHeuristicRandRound(const CbcHeuristicRandRound & rhs)
66:
67  CbcHeuristic(rhs)
68{
69}
70
71// Assignment operator
72CbcHeuristicRandRound & 
73CbcHeuristicRandRound::operator=( const CbcHeuristicRandRound& rhs)
74{
75  if (this!=&rhs) {
76    CbcHeuristic::operator=(rhs);
77  }
78  return *this;
79}
80
81// Resets stuff if model changes
82void 
83CbcHeuristicRandRound::resetModel(CbcModel * model)
84{
85  CbcHeuristic::resetModel(model);
86}
87
88/*
89  Randomized Rounding Heuristic
90  Returns 1 if solution, 0 if not
91*/
92int
93CbcHeuristicRandRound::solution(double & solutionValue,
94                         double * betterSolution)
95{
96  // rlh: Todo: Memory Cleanup
97
98  //  std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;
99
100  setWhen(1);  // setWhen(1) didn't have the effect I expected (e.g., run once).
101             
102  // Run only once.
103  //
104  //    See if at root node
105  bool atRoot = model_->getNodeCount()==0;
106  int passNumber = model_->getCurrentPassNumber();
107  //    Just do once
108  if (!atRoot||passNumber!=1){
109    // std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
110    return 0;
111  }
112   
113  std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;
114  typedef struct {
115    int numberSolutions;
116    int maximumSolutions;
117    int numberColumns;
118    double ** solution;
119    int * numberUnsatisfied;
120  } clpSolution;
121
122  double start = CoinCpuTime();
123  numCouldRun_++; //
124                  // Todo: Ask JJHF what "number of times
125                  // the heuristic could run" means.
126
127  OsiSolverInterface * solver= model_->solver()->clone();
128  OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *> (solver);
129  assert (clpSolver);
130  ClpSimplex * simplex = clpSolver->getModelPtr();
131
132  // Initialize the structure holding the solutions for the Simplex iterations
133  clpSolution solutions;
134  // Set typeStruct field of ClpTrustedData struct to 1 to indicate
135  // desired behavior for  RandRound heuristic (which is what?)
136  ClpTrustedData trustedSolutions;
137  trustedSolutions.typeStruct = 1;
138  trustedSolutions.data = &solutions;
139  solutions.numberSolutions=0;
140  solutions.maximumSolutions=0;
141  solutions.numberColumns=simplex->numberColumns();
142  solutions.solution=NULL;
143  solutions.numberUnsatisfied=NULL;
144  simplex->setTrustedUserPointer(&trustedSolutions);
145
146  // Solve from all slack to get some points
147  simplex->allSlackBasis();
148 
149  // Calling primal() invalidates pointers to some rim vectors,
150  // like...row sense (!)
151  simplex->primal();
152
153  // 1. Okay - so a workaround would be to copy the data I want BEFORE
154  // calling primal. 
155  // 2. Another approach is to ask the simplex solvers NOT to mess up my
156  // rims.
157  // 3. See freeCachedResults() for what is getting
158  // deleted. Everything else points into the structure.
159  // ...or use collower and colupper rather than rowsense.
160  // ..store address of where one of these
161
162  // Store the basic problem information
163  // -Get the number of columns, rows and rhs vector
164  int numCols = clpSolver->getNumCols();
165  int numRows = clpSolver->getNumRows();
166
167  // Find the integer variables (use columnType(?))
168  // One if not continuous, that is binary or general integer)
169  // columnType() = 0 continuous
170  //              = 1 binary
171  //              = 2 general integer
172  bool * varClassInt = new bool[numCols];
173  const char* columnType = clpSolver->columnType();
174  int numGenInt = 0;   
175  for(int i=0; i<numCols; i++)
176    {
177      if(clpSolver->isContinuous(i))
178        varClassInt[i] = 0;
179      else
180        varClassInt[i] = 1;
181      if (columnType[i]==2) numGenInt++;
182    }
183
184  // Heuristic is for problems with general integer variables.
185  // If there are none, quit.
186  if (numGenInt++<1){
187    delete [] varClassInt ;
188    std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
189    return 0; 
190 }
191
192
193  // -Get the rows sense
194  const char * rowSense;
195  rowSense = clpSolver->getRowSense();
196
197  // -Get the objective coefficients
198  double *originalObjCoeff = CoinCopyOfArray(clpSolver->getObjCoefficients(), numCols);
199
200  // -Get the matrix of the problem
201  // rlh: look at using sparse representation
202 double ** matrix = new double * [numRows];
203  for(int i = 0; i < numRows; i++)
204    {
205      matrix[i] = new double[numCols];
206      for(int j = 0; j < numCols; j++)
207        matrix[i][j] = 0;
208    }
209 
210  const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
211  const double * matrixElements = matrixByRow->getElements();
212  const int * matrixIndices = matrixByRow->getIndices();
213  const int * matrixStarts = matrixByRow->getVectorStarts();
214  for(int j = 0; j < numRows; j++)
215    {
216      for(int i = matrixStarts[j]; i < matrixStarts[j+1]; i++)
217        {
218          matrix[j][matrixIndices[i]] = matrixElements[i];
219        }
220    }
221 
222  double * newObj = new double [numCols];
223  srand ( time(NULL) +1);
224  int randNum;
225
226  // Shuffle the rows:
227  // Put the rows in a random order
228  // so that the optimal solution is a different corner point than the
229  // starting point.
230  int * index = new int [numRows];
231  for(int i=0; i<numRows; i++)
232    index[i]=i;
233  for(int i=0; i<numRows; i++)
234    {
235      int temp = index[i];
236      int randNumTemp = i+intRand(numRows-i);
237      index[i]=index[randNumTemp];
238      index[randNumTemp]=temp;
239    }
240 
241  // Start finding corner points by iteratively doing the following:
242  // - contruct a randomly tilted objective
243  // - solve
244  for(int i=0; i<numRows; i++)
245    {
246      // TODO: that 10,000 could be a param in the member data
247      if(solutions.numberSolutions  > 10000)
248        break;
249      randNum = intRand(2);
250      for(int j=0; j<numCols; j++)
251        {
252          // for row i and column j vary the coefficient "a bit"
253          if(randNum == 1)             
254            // if the element is zero, then set the new obj
255            // coefficient to 0.1 (i.e., round up)
256            if(fabs(matrix[index[i]][j]) < 1e-6)
257              newObj[j] = 0.1;
258            else
259              // if the element is nonzero, then increase the new obj
260              // coefficient "a bit"
261              newObj[j] = matrix[index[i]][j] * 1.1;
262          else
263            // if randnum is 2, then
264            // if the element is zero, then set the new obj coeffient
265            // to NEGATIVE 0.1 (i.e., round down)
266            if(fabs(matrix[index[i]][j]) < 1e-6)
267              newObj[j] = -0.1;
268            else
269              // if the element is nonzero, then DEcrease the new obj coeffienct "a bit"
270              newObj[j] = matrix[index[i]][j] * 0.9;
271        }
272      // Use the new "tilted" objective
273      clpSolver->setObjective(newObj);
274
275      // Based on the row sense, we decide whether to max or min
276      if(rowSense[i] == 'L')
277        clpSolver->setObjSense(-1);
278      else
279        clpSolver->setObjSense(1);
280
281      // Solve with primal simplex
282      simplex->primal(1);
283      // rlh+ll: This was the original code. But we already have the
284      // model pointer (it's in simplex). And, calling getModelPtr()
285      // invalidates the cached data in the OsiClpSolverInterface
286      // object, which means our precious rowsens is lost. So let's
287      // not use the line below...
288      /******* clpSolver->getModelPtr()->primal(1); */
289      printf("---------------------------------------------------------------- %d\n", i);
290    }
291  // Iteratively do this process until...
292  // either you reach the max number of corner points (aka 10K)
293  // or all the rows have been used as an objective.
294 
295  // Look at solutions
296  int numberSolutions = solutions.numberSolutions;
297  //const char * integerInfo = simplex->integerInformation();
298  //const double * columnLower = simplex->columnLower();
299  //const double * columnUpper = simplex->columnUpper();
300  printf("there are %d solutions\n",numberSolutions);
301
302  // Up to here we have all the corner points
303  // Now we need to do the random walks and roundings
304 
305  double ** cornerPoints = new double * [numberSolutions];
306  for(int j=0;j<numberSolutions;j++)
307    cornerPoints[j] = solutions.solution[j];
308 
309  bool feasibility = 1;
310  // rlh: use some COIN max instead of 1e30 (?)
311  double bestObj = 1e30;
312  std::vector< std::vector <double> > feasibles;
313  int numFeasibles = 0;
314 
315  // Check the feasibility of the corner points
316  int numCornerPoints = numberSolutions;
317
318  const double * rhs = clpSolver->getRightHandSide();
319  // rlh: row sense hasn't changed. why a fresh copy?
320  // Delete next line.
321  rowSense = clpSolver->getRowSense();
322 
323  for(int i=0; i<numCornerPoints; i++)
324    {
325      //get the objective value for this this point
326      double objValue = 0;
327      for(int k=0; k<numCols; k++)
328        objValue += cornerPoints[i][k] * originalObjCoeff[k];
329     
330      if(objValue < bestObj)
331        {
332          // check integer feasibility
333          feasibility = 1;
334          for(int j=0; j<numCols; j++)
335            {
336              if(varClassInt[j])
337                {
338                  double closest=floor(cornerPoints[i][j]+0.5);
339                  if(fabs(cornerPoints[i][j] - closest)>1e-6)
340                    {
341                      feasibility = 0;
342                      break;
343                    }
344                }
345            }
346          // check all constraints satisfied
347          if (feasibility)
348              { 
349                for(int irow = 0; irow < numRows; irow++)
350                  {
351                    double lhs = 0;
352                    for(int j = 0; j <numCols; j++)
353                      {
354                        lhs += matrix[irow][j] * cornerPoints[i][j];
355                      }
356                    if(rowSense[irow] == 'L' && lhs > rhs[irow] + 1e-6)
357                      {
358                        feasibility = 0;
359                        break;
360                      }
361                    if(rowSense[irow] == 'G' && lhs < rhs[irow] - 1e-6)
362                      {
363                        feasibility = 0;
364                        break;
365                      }
366                    if(rowSense[irow] == 'E' && (lhs - rhs[irow] > 1e-6 || lhs - rhs[irow] < -1e-6)) 
367                      {
368                        feasibility = 0;
369                        break;
370                      }
371                  }
372              }
373         
374          if(feasibility)
375            {           
376              numFeasibles++;
377                feasibles.push_back(std::vector <double> (numCols));
378                for(int k=0; k<numCols; k++)
379                  feasibles[numFeasibles-1][k] = cornerPoints[i][k];
380                printf("obj: %f\n", objValue);
381                if(objValue < bestObj)
382                  bestObj = objValue;
383            }
384        }
385    }   
386  int numFeasibleCorners;
387  numFeasibleCorners = numFeasibles;
388  //find the center of gravity of the corner points as the first random point
389  double * rp = new double[numCols];
390  for(int i=0; i<numCols; i++)
391    {
392      rp[i] = 0;
393      for(int j=0; j < numCornerPoints; j++)
394        {
395            rp[i] += cornerPoints[j][i];
396        }
397      rp[i] = rp[i]/numCornerPoints;
398    }
399 
400  //-------------------------------------------
401  //main loop:
402  // -generate the next random point
403  // -round the random point
404  // -check the feasibility of the random point
405  //-------------------------------------------
406 
407  srand ( time(NULL) +1);
408  int numRandomPoints = 0;
409  while(numRandomPoints < 50000)
410    {
411      numRandomPoints++;
412      //generate the next random point
413      int randomIndex = intRand(numCornerPoints);
414      double random = CoinDrand48();
415      for(int i=0; i<numCols; i++)
416        {
417          rp[i] = (random * (cornerPoints[randomIndex][i] - rp[i])) + rp[i];
418        }
419     
420      //CRISP ROUNDING
421      //round the random point just generated
422      double * roundRp = new double[numCols];
423      for(int i=0; i<numCols; i++)
424        {
425          roundRp[i] = rp[i];
426          if(varClassInt[i])
427            {
428              if(rp[i]>=0)
429                {                               
430                  if(fmod(rp[i],1) > 0.5)
431                    roundRp[i] = floor(rp[i]) + 1;
432                  else
433                    roundRp[i] = floor(rp[i]);
434                }
435              else
436                {
437                  if(fabs(fmod(rp[i],1)) > 0.5)
438                    roundRp[i] = floor(rp[i]);
439                  else
440                    roundRp[i] = floor(rp[i])+1;
441                 
442                }
443            }
444        }
445     
446     
447      //SOFT ROUNDING
448      // Look at original files for the "how to" on soft rounding;
449      // Soft rounding omitted here. 
450           
451      //Check the feasibility of the rounded random point
452      // -Check the feasibility
453      // -Get the rows sense
454      rowSense = clpSolver->getRowSense();
455      rhs = clpSolver->getRightHandSide();
456     
457      //get the objective value for this feasible point
458      double objValue = 0;
459      for(int i=0; i<numCols; i++)
460        objValue += roundRp[i] * originalObjCoeff[i];
461     
462      if(objValue<bestObj)
463        {
464          feasibility = 1;
465          for(int i = 0; i < numRows; i++)
466            {
467              double lhs = 0;
468              for(int j = 0; j <numCols; j++)
469                {
470                  lhs += matrix[i][j] * roundRp[j];
471                }
472              if(rowSense[i] == 'L' && lhs > rhs[i] + 1e-6)
473                {
474                  feasibility = 0;
475                  break;
476                }
477              if(rowSense[i] == 'G' && lhs < rhs[i] - 1e-6)
478                {
479                  feasibility = 0;
480                  break;
481                }
482              if(rowSense[i] == 'E' && (lhs - rhs[i] > 1e-6 || lhs - rhs[i] < -1e-6)) 
483                {
484                  feasibility = 0;
485                  break;
486                }
487            }
488          if(feasibility)
489            {
490              printf("Feasible Found.\n");
491              printf("%.2f\n", CoinCpuTime()-start);
492              numFeasibles++;
493              feasibles.push_back(std::vector <double> (numCols));
494              for(int i=0; i<numCols; i++)
495                feasibles[numFeasibles-1][i] = roundRp[i];
496              printf("obj: %f\n", objValue);
497              if(objValue < bestObj)
498                bestObj = objValue;
499            }
500        }
501      delete [] roundRp;
502    }
503  printf("Number of Feasible Corners: %d\n", numFeasibleCorners);
504  printf("Number of Feasibles Found: %d\n", numFeasibles);
505  if(numFeasibles > 0)
506    printf("Best Objective: %f\n", bestObj);
507  printf("time: %.2f\n",CoinCpuTime()-start);
508 
509  if (numFeasibles == 0) {
510    // cleanup
511    delete [] varClassInt;
512    for (int i=0; i<numRows; i++)
513      delete matrix[i];
514    delete [] matrix; 
515    delete [] newObj;
516    delete [] index;
517    for (int i=0; i<numberSolutions; i++)
518      delete cornerPoints[i];
519    delete [] cornerPoints; 
520    delete [] rp;
521    return 0; 
522  }
523 
524  // We found something better
525  solutionValue = bestObj;
526  for (int k = 0; k<numCols; k++){
527    betterSolution[k] =  feasibles[numFeasibles-1][k];
528  } 
529  delete [] varClassInt;
530  for (int i=0; i<numRows; i++)
531    delete matrix[i];
532  delete [] matrix; 
533  delete [] newObj;
534  delete [] index;
535  for (int i=0; i<numberSolutions; i++)
536      delete cornerPoints[i];
537    delete [] cornerPoints; 
538  delete [] rp;
539  std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
540  return 1;
541
542}
543// update model
544void CbcHeuristicRandRound::setModel(CbcModel * model)
545{
546  CbcHeuristic::setModel(model);
547}
548
549 
Note: See TracBrowser for help on using the repository browser.