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