source: stable/2.4/Cbc/src/CbcHeuristicRandRound.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

File size: 15.1 KB
Line 
1/* $Id: CbcHeuristicRandRound.cpp 1200 2009-07-25 08:44:13Z forrest $ */
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 static_cast<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  double primalTolerance ;
129  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
130  OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *> (solver);
131  assert (clpSolver);
132  ClpSimplex * simplex = clpSolver->getModelPtr();
133
134  // Initialize the structure holding the solutions for the Simplex iterations
135  clpSolution solutions;
136  // Set typeStruct field of ClpTrustedData struct to 1 to indicate
137  // desired behavior for  RandRound heuristic (which is what?)
138  ClpTrustedData trustedSolutions;
139  trustedSolutions.typeStruct = 1;
140  trustedSolutions.data = &solutions;
141  solutions.numberSolutions=0;
142  solutions.maximumSolutions=0;
143  solutions.numberColumns=simplex->numberColumns();
144  solutions.solution=NULL;
145  solutions.numberUnsatisfied=NULL;
146  simplex->setTrustedUserPointer(&trustedSolutions);
147
148  // Solve from all slack to get some points
149  simplex->allSlackBasis();
150 
151  // Calling primal() invalidates pointers to some rim vectors,
152  // like...row sense (!)
153  simplex->primal();
154
155  // 1. Okay - so a workaround would be to copy the data I want BEFORE
156  // calling primal. 
157  // 2. Another approach is to ask the simplex solvers NOT to mess up my
158  // rims.
159  // 3. See freeCachedResults() for what is getting
160  // deleted. Everything else points into the structure.
161  // ...or use collower and colupper rather than rowsense.
162  // ..store address of where one of these
163
164  // Store the basic problem information
165  // -Get the number of columns, rows and rhs vector
166  int numCols = clpSolver->getNumCols();
167  int numRows = clpSolver->getNumRows();
168
169  // Find the integer variables (use columnType(?))
170  // One if not continuous, that is binary or general integer)
171  // columnType() = 0 continuous
172  //              = 1 binary
173  //              = 2 general integer
174  bool * varClassInt = new bool[numCols];
175  const char* columnType = clpSolver->columnType();
176  int numGenInt = 0;   
177  for(int i=0; i<numCols; i++)
178    {
179      if(clpSolver->isContinuous(i))
180        varClassInt[i] = 0;
181      else
182        varClassInt[i] = 1;
183      if (columnType[i]==2) numGenInt++;
184    }
185
186  // Heuristic is for problems with general integer variables.
187  // If there are none, quit.
188  if (numGenInt++<1){
189    delete [] varClassInt ;
190    std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
191    return 0; 
192 }
193
194
195  // -Get the rows sense
196  const char * rowSense;
197  rowSense = clpSolver->getRowSense();
198
199  // -Get the objective coefficients
200  double *originalObjCoeff = CoinCopyOfArray(clpSolver->getObjCoefficients(), numCols);
201
202  // -Get the matrix of the problem
203  // rlh: look at using sparse representation
204 double ** matrix = new double * [numRows];
205  for(int i = 0; i < numRows; i++)
206    {
207      matrix[i] = new double[numCols];
208      for(int j = 0; j < numCols; j++)
209        matrix[i][j] = 0;
210    }
211 
212  const CoinPackedMatrix* matrixByRow = clpSolver->getMatrixByRow();
213  const double * matrixElements = matrixByRow->getElements();
214  const int * matrixIndices = matrixByRow->getIndices();
215  const int * matrixStarts = matrixByRow->getVectorStarts();
216  for(int j = 0; j < numRows; j++)
217    {
218      for(int i = matrixStarts[j]; i < matrixStarts[j+1]; i++)
219        {
220          matrix[j][matrixIndices[i]] = matrixElements[i];
221        }
222    }
223 
224  double * newObj = new double [numCols];
225  srand ( time(NULL) +1);
226  int randNum;
227
228  // Shuffle the rows:
229  // Put the rows in a random order
230  // so that the optimal solution is a different corner point than the
231  // starting point.
232  int * index = new int [numRows];
233  for(int i=0; i<numRows; i++)
234    index[i]=i;
235  for(int i=0; i<numRows; i++)
236    {
237      int temp = index[i];
238      int randNumTemp = i+intRand(numRows-i);
239      index[i]=index[randNumTemp];
240      index[randNumTemp]=temp;
241    }
242 
243  // Start finding corner points by iteratively doing the following:
244  // - contruct a randomly tilted objective
245  // - solve
246  for(int i=0; i<numRows; i++)
247    {
248      // TODO: that 10,000 could be a param in the member data
249      if(solutions.numberSolutions  > 10000)
250        break;
251      randNum = intRand(2);
252      for(int j=0; j<numCols; j++)
253        {
254          // for row i and column j vary the coefficient "a bit"
255          if(randNum == 1)             
256            // if the element is zero, then set the new obj
257            // coefficient to 0.1 (i.e., round up)
258            if(fabs(matrix[index[i]][j]) < primalTolerance)
259              newObj[j] = 0.1;
260            else
261              // if the element is nonzero, then increase the new obj
262              // coefficient "a bit"
263              newObj[j] = matrix[index[i]][j] * 1.1;
264          else
265            // if randnum is 2, then
266            // if the element is zero, then set the new obj coeffient
267            // to NEGATIVE 0.1 (i.e., round down)
268            if(fabs(matrix[index[i]][j]) < primalTolerance)
269              newObj[j] = -0.1;
270            else
271              // if the element is nonzero, then DEcrease the new obj coeffienct "a bit"
272              newObj[j] = matrix[index[i]][j] * 0.9;
273        }
274      // Use the new "tilted" objective
275      clpSolver->setObjective(newObj);
276
277      // Based on the row sense, we decide whether to max or min
278      if(rowSense[i] == 'L')
279        clpSolver->setObjSense(-1);
280      else
281        clpSolver->setObjSense(1);
282
283      // Solve with primal simplex
284      simplex->primal(1);
285      // rlh+ll: This was the original code. But we already have the
286      // model pointer (it's in simplex). And, calling getModelPtr()
287      // invalidates the cached data in the OsiClpSolverInterface
288      // object, which means our precious rowsens is lost. So let's
289      // not use the line below...
290      /******* clpSolver->getModelPtr()->primal(1); */
291      printf("---------------------------------------------------------------- %d\n", i);
292    }
293  // Iteratively do this process until...
294  // either you reach the max number of corner points (aka 10K)
295  // or all the rows have been used as an objective.
296 
297  // Look at solutions
298  int numberSolutions = solutions.numberSolutions;
299  //const char * integerInfo = simplex->integerInformation();
300  //const double * columnLower = simplex->columnLower();
301  //const double * columnUpper = simplex->columnUpper();
302  printf("there are %d solutions\n",numberSolutions);
303
304  // Up to here we have all the corner points
305  // Now we need to do the random walks and roundings
306 
307  double ** cornerPoints = new double * [numberSolutions];
308  for(int j=0;j<numberSolutions;j++)
309    cornerPoints[j] = solutions.solution[j];
310 
311  bool feasibility = 1;
312  // rlh: use some COIN max instead of 1e30 (?)
313  double bestObj = 1e30;
314  std::vector< std::vector <double> > feasibles;
315  int numFeasibles = 0;
316 
317  // Check the feasibility of the corner points
318  int numCornerPoints = numberSolutions;
319
320  const double * rhs = clpSolver->getRightHandSide();
321  // rlh: row sense hasn't changed. why a fresh copy?
322  // Delete next line.
323  rowSense = clpSolver->getRowSense();
324 
325  for(int i=0; i<numCornerPoints; i++)
326    {
327      //get the objective value for this this point
328      double objValue = 0;
329      for(int k=0; k<numCols; k++)
330        objValue += cornerPoints[i][k] * originalObjCoeff[k];
331     
332      if(objValue < bestObj)
333        {
334          // check integer feasibility
335          feasibility = 1;
336          for(int j=0; j<numCols; j++)
337            {
338              if(varClassInt[j])
339                {
340                  double closest=floor(cornerPoints[i][j]+0.5);
341                  if(fabs(cornerPoints[i][j] - closest)>primalTolerance)
342                    {
343                      feasibility = 0;
344                      break;
345                    }
346                }
347            }
348          // check all constraints satisfied
349          if (feasibility)
350              { 
351                for(int irow = 0; irow < numRows; irow++)
352                  {
353                    double lhs = 0;
354                    for(int j = 0; j <numCols; j++)
355                      {
356                        lhs += matrix[irow][j] * cornerPoints[i][j];
357                      }
358                    if(rowSense[irow] == 'L' && lhs > rhs[irow] + primalTolerance)
359                      {
360                        feasibility = 0;
361                        break;
362                      }
363                    if(rowSense[irow] == 'G' && lhs < rhs[irow] - primalTolerance)
364                      {
365                        feasibility = 0;
366                        break;
367                      }
368                    if(rowSense[irow] == 'E' && (lhs - rhs[irow] > primalTolerance || lhs - rhs[irow] < -primalTolerance)) 
369                      {
370                        feasibility = 0;
371                        break;
372                      }
373                  }
374              }
375         
376          if(feasibility)
377            {           
378              numFeasibles++;
379                feasibles.push_back(std::vector <double> (numCols));
380                for(int k=0; k<numCols; k++)
381                  feasibles[numFeasibles-1][k] = cornerPoints[i][k];
382                printf("obj: %f\n", objValue);
383                if(objValue < bestObj)
384                  bestObj = objValue;
385            }
386        }
387    }   
388  int numFeasibleCorners;
389  numFeasibleCorners = numFeasibles;
390  //find the center of gravity of the corner points as the first random point
391  double * rp = new double[numCols];
392  for(int i=0; i<numCols; i++)
393    {
394      rp[i] = 0;
395      for(int j=0; j < numCornerPoints; j++)
396        {
397            rp[i] += cornerPoints[j][i];
398        }
399      rp[i] = rp[i]/numCornerPoints;
400    }
401 
402  //-------------------------------------------
403  //main loop:
404  // -generate the next random point
405  // -round the random point
406  // -check the feasibility of the random point
407  //-------------------------------------------
408 
409  srand ( time(NULL) +1);
410  int numRandomPoints = 0;
411  while(numRandomPoints < 50000)
412    {
413      numRandomPoints++;
414      //generate the next random point
415      int randomIndex = intRand(numCornerPoints);
416      double random = CoinDrand48();
417      for(int i=0; i<numCols; i++)
418        {
419          rp[i] = (random * (cornerPoints[randomIndex][i] - rp[i])) + rp[i];
420        }
421     
422      //CRISP ROUNDING
423      //round the random point just generated
424      double * roundRp = new double[numCols];
425      for(int i=0; i<numCols; i++)
426        {
427          roundRp[i] = rp[i];
428          if(varClassInt[i])
429            {
430              if(rp[i]>=0)
431                {                               
432                  if(fmod(rp[i],1) > 0.5)
433                    roundRp[i] = floor(rp[i]) + 1;
434                  else
435                    roundRp[i] = floor(rp[i]);
436                }
437              else
438                {
439                  if(fabs(fmod(rp[i],1)) > 0.5)
440                    roundRp[i] = floor(rp[i]);
441                  else
442                    roundRp[i] = floor(rp[i])+1;
443                 
444                }
445            }
446        }
447     
448     
449      //SOFT ROUNDING
450      // Look at original files for the "how to" on soft rounding;
451      // Soft rounding omitted here. 
452           
453      //Check the feasibility of the rounded random point
454      // -Check the feasibility
455      // -Get the rows sense
456      rowSense = clpSolver->getRowSense();
457      rhs = clpSolver->getRightHandSide();
458     
459      //get the objective value for this feasible point
460      double objValue = 0;
461      for(int i=0; i<numCols; i++)
462        objValue += roundRp[i] * originalObjCoeff[i];
463     
464      if(objValue<bestObj)
465        {
466          feasibility = 1;
467          for(int i = 0; i < numRows; i++)
468            {
469              double lhs = 0;
470              for(int j = 0; j <numCols; j++)
471                {
472                  lhs += matrix[i][j] * roundRp[j];
473                }
474              if(rowSense[i] == 'L' && lhs > rhs[i] + primalTolerance)
475                {
476                  feasibility = 0;
477                  break;
478                }
479              if(rowSense[i] == 'G' && lhs < rhs[i] - primalTolerance)
480                {
481                  feasibility = 0;
482                  break;
483                }
484              if(rowSense[i] == 'E' && (lhs - rhs[i] > primalTolerance || lhs - rhs[i] < -primalTolerance)) 
485                {
486                  feasibility = 0;
487                  break;
488                }
489            }
490          if(feasibility)
491            {
492              printf("Feasible Found.\n");
493              printf("%.2f\n", CoinCpuTime()-start);
494              numFeasibles++;
495              feasibles.push_back(std::vector <double> (numCols));
496              for(int i=0; i<numCols; i++)
497                feasibles[numFeasibles-1][i] = roundRp[i];
498              printf("obj: %f\n", objValue);
499              if(objValue < bestObj)
500                bestObj = objValue;
501            }
502        }
503      delete [] roundRp;
504    }
505  printf("Number of Feasible Corners: %d\n", numFeasibleCorners);
506  printf("Number of Feasibles Found: %d\n", numFeasibles);
507  if(numFeasibles > 0)
508    printf("Best Objective: %f\n", bestObj);
509  printf("time: %.2f\n",CoinCpuTime()-start);
510 
511  if (numFeasibles == 0) {
512    // cleanup
513    delete [] varClassInt;
514    for (int i=0; i<numRows; i++)
515      delete matrix[i];
516    delete [] matrix; 
517    delete [] newObj;
518    delete [] index;
519    for (int i=0; i<numberSolutions; i++)
520      delete cornerPoints[i];
521    delete [] cornerPoints; 
522    delete [] rp;
523    return 0; 
524  }
525 
526  // We found something better
527  solutionValue = bestObj;
528  for (int k = 0; k<numCols; k++){
529    betterSolution[k] =  feasibles[numFeasibles-1][k];
530  } 
531  delete [] varClassInt;
532  for (int i=0; i<numRows; i++)
533    delete matrix[i];
534  delete [] matrix; 
535  delete [] newObj;
536  delete [] index;
537  for (int i=0; i<numberSolutions; i++)
538      delete cornerPoints[i];
539    delete [] cornerPoints; 
540  delete [] rp;
541  std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
542  return 1;
543
544}
545// update model
546void CbcHeuristicRandRound::setModel(CbcModel * model)
547{
548  CbcHeuristic::setModel(model);
549}
550
551 
Note: See TracBrowser for help on using the repository browser.