source: trunk/Cbc/src/CbcHeuristicPivotAndFix.cpp @ 2464

Last change on this file since 2464 was 2464, checked in by unxusr, 6 months ago

.clang-format with proposal for formatting code

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