source: trunk/Cbc/src/CbcHeuristicRandRound.cpp

Last change on this file was 2467, checked in by unxusr, 5 months ago

spaces after angles

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