source: trunk/Cbc/examples/allCuts.cpp @ 1565

Last change on this file since 1565 was 1468, checked in by stefan, 9 years ago

do not require CbcConfig?.h in example to decide whether sample or miplib3 is present - do this in makefile

  • Property svn:executable set to *
File size: 16.8 KB
Line 
1// Copyright (C) 2006, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include <cassert>
9#include <iomanip>
10
11
12// For Branch and bound
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcBranchUser.hpp"
16#include "CbcCompareUser.hpp"
17#include "CbcCutGenerator.hpp"
18#include "CbcStrategy.hpp"
19#include "CbcHeuristicLocal.hpp"
20#include "OsiClpSolverInterface.hpp"
21
22// Cuts
23
24#include "CglGomory.hpp"
25#include "CglProbing.hpp"
26#include "CglKnapsackCover.hpp"
27#include "CglRedSplit.hpp"
28#include "CglClique.hpp"
29#include "CglFlowCover.hpp"
30#include "CglMixedIntegerRounding2.hpp"
31// Preprocessing
32#include "CglPreProcess.hpp"
33
34// For saying about solution validity
35#include "OsiAuxInfo.hpp"
36
37// Heuristics (But any will have to be special)
38
39#include "CbcHeuristic.hpp"
40
41#include  "CoinTime.hpp"
42// Need stored cuts
43
44#include "CglStored.hpp"
45
46/** Stored Cut Generator Class */
47class CglStoredUser : public CglStored {
48 
49public:
50   
51 
52  /**@name Generate Cuts */
53  //@{
54  /** Generate Mixed Integer Stored cuts for the model of the
55      solver interface, si.
56
57      Insert the generated cuts into OsiCut, cs.
58
59      This generator just looks at previously stored cuts
60      and inserts any that are violated by enough
61  */
62  virtual void generateCuts( const OsiSolverInterface & si, OsiCuts & cs,
63                             const CglTreeInfo info = CglTreeInfo()) const;
64  //@}
65
66  /**@name Cut stuff */
67  //@{
68  OsiRowCut * mutableRowCutPointer(int index)
69  { return cuts_.rowCutPtr(index);}
70  //@}
71
72  /**@name Constructors and destructors */
73  //@{
74  /// Default constructor
75  CglStoredUser ();
76 
77  /// Copy constructor
78  CglStoredUser (const CglStoredUser & rhs);
79
80  /// Clone
81  virtual CglCutGenerator * clone() const;
82
83  /// Assignment operator
84  CglStoredUser &
85    operator=(const CglStoredUser& rhs);
86 
87  /// Destructor
88  virtual
89    ~CglStoredUser ();
90  //@}
91     
92protected:
93 
94 // Protected member methods
95
96  // Protected member data
97
98  /**@name Protected member data */
99  //@{
100  /** Don't add any more cuts after this number passes (per node)
101      unless looks integer feasible
102  */
103  int numberPasses_;
104  //@}
105};
106//-------------------------------------------------------------------
107// Generate Stored cuts
108//-------------------------------------------------------------------
109void 
110CglStoredUser::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
111                             const CglTreeInfo info) const
112{
113  // Get basic problem information
114  const double * solution = si.getColSolution();
115  if (info.inTree&&info.pass>numberPasses_) {
116    // only continue if integer feasible
117    int numberColumns=si.getNumCols(); 
118    int i;
119    const double * colUpper = si.getColUpper();
120    const double * colLower = si.getColLower();
121    int numberAway=0;
122    for (i=0;i<numberColumns;i++) {
123      double value = solution[i];
124      // In case slightly away from bounds
125      value = CoinMax(colLower[i],value);
126      value = CoinMin(colUpper[i],value);
127      if (si.isInteger(i)&&fabs(value-fabs(value+0.5))>1.0e-5) 
128        numberAway++;
129    }
130    if (numberAway)
131      return; // let code branch
132  }
133  int numberRowCuts = cuts_.sizeRowCuts();
134  for (int i=0;i<numberRowCuts;i++) {
135    const OsiRowCut * rowCutPointer = cuts_.rowCutPtr(i);
136    double violation = rowCutPointer->violated(solution);
137    if (violation>=requiredViolation_)
138      cs.insert(*rowCutPointer);
139  }
140}
141
142//-------------------------------------------------------------------
143// Default Constructor
144//-------------------------------------------------------------------
145CglStoredUser::CglStoredUser ()
146:
147  CglStored(),
148  numberPasses_(5)
149{
150}
151
152//-------------------------------------------------------------------
153// Copy constructor
154//-------------------------------------------------------------------
155CglStoredUser::CglStoredUser (const CglStoredUser & source) :
156  CglStored(source),
157  numberPasses_(source.numberPasses_)
158{ 
159}
160
161//-------------------------------------------------------------------
162// Clone
163//-------------------------------------------------------------------
164CglCutGenerator *
165CglStoredUser::clone() const
166{
167  return new CglStoredUser(*this);
168}
169
170//-------------------------------------------------------------------
171// Destructor
172//-------------------------------------------------------------------
173CglStoredUser::~CglStoredUser ()
174{
175}
176
177//----------------------------------------------------------------
178// Assignment operator
179//-------------------------------------------------------------------
180CglStoredUser &
181CglStoredUser::operator=(const CglStoredUser& rhs)
182{
183  if (this != &rhs) {
184    CglStored::operator=(rhs);
185    numberPasses_=rhs.numberPasses_;
186  }
187  return *this;
188}
189// Class to disallow strong branching solutions
190#include "CbcFeasibilityBase.hpp"
191class CbcFeasibilityNoStrong : public CbcFeasibilityBase{
192public:
193  // Default Constructor
194  CbcFeasibilityNoStrong () {};
195
196  virtual ~CbcFeasibilityNoStrong() {};
197  // Copy constructor
198  CbcFeasibilityNoStrong ( const CbcFeasibilityNoStrong &rhs) {};
199   
200  // Assignment operator
201  CbcFeasibilityNoStrong & operator=( const CbcFeasibilityNoStrong& rhs)
202  { return * this;};
203
204  /// Clone
205  virtual CbcFeasibilityBase * clone() const
206  { return new CbcFeasibilityNoStrong();};
207
208  /**
209     On input mode:
210     0 - called after a solve but before any cuts
211     -1 - called after strong branching
212     Returns :
213     0 - no opinion
214     -1 pretend infeasible
215     1 pretend integer solution
216  */
217  virtual int feasible(CbcModel * model, int mode) 
218  {return mode;};
219};
220
221
222//#############################################################################
223
224
225/************************************************************************
226
227This main program reads in an integer model from an mps file.
228
229It makes L or G rows into cuts
230
231************************************************************************/
232
233int main (int argc, const char *argv[])
234{
235
236  // Define your favorite OsiSolver
237 
238  OsiClpSolverInterface solver1;
239
240  // Read in model using argv[1]
241  // and assert that it is a clean model
242  std::string mpsFileName;
243#if defined(SAMPLEDIR)
244  mpsFileName = SAMPLEDIR "/p0033.mps";
245#else
246  if (argc < 2) {
247    fprintf(stderr, "Do not know where to find sample MPS files.\n");
248    exit(1);
249  }
250#endif
251  if (argc>=2) mpsFileName = argv[1];
252  int numMpsReadErrors = solver1.readMps(mpsFileName.c_str(),"");
253  assert(numMpsReadErrors==0);
254  double time1 = CoinCpuTime();
255  OsiClpSolverInterface solverSave = solver1;
256
257  /* Options are:
258     preprocess to do preprocessing
259     time in minutes
260     if 2 parameters and numeric taken as time
261  */
262  bool preProcess=false;
263  double minutes=-1.0;
264  int nGoodParam=0;
265  for (int iParam=2; iParam<argc;iParam++) {
266    if (!strcmp(argv[iParam],"preprocess")) {
267      preProcess=true;
268      nGoodParam++;
269    } else if (!strcmp(argv[iParam],"time")) {
270      if (iParam+1<argc&&isdigit(argv[iParam+1][0])) {
271        minutes=atof(argv[iParam+1]);
272        if (minutes>=0.0) {
273          nGoodParam+=2;
274          iParam++; // skip time
275        }
276      }
277    }
278  }
279  if (nGoodParam==0&&argc==3&&isdigit(argv[2][0])) {
280    // If time is given then stop after that number of minutes
281    minutes = atof(argv[2]);
282    if (minutes>=0.0) 
283      nGoodParam=1;
284  }
285  if (nGoodParam!=argc-2&&argc>=2) {
286    printf("Usage <file> [preprocess] [time <minutes>] or <file> <minutes>\n");
287    exit(1);
288  }
289  // Reduce printout
290  solver1.setHintParam(OsiDoReducePrint,true,OsiHintTry);
291  // See if we want preprocessing
292  OsiSolverInterface * solver2=&solver1;
293  CglPreProcess process;
294  // Never do preprocessing until dual tests out as can fix incorrectly
295  preProcess=false;
296  if (preProcess) {
297    /* Do not try and produce equality cliques and
298       do up to 5 passes */
299    solver2 = process.preProcess(solver1,false,5);
300    if (!solver2) {
301      printf("Pre-processing says infeasible\n");
302      exit(2);
303    }
304    solver2->resolve();
305  }
306  // Turn L rows into cuts
307  CglStoredUser stored;
308 {
309  int numberRows = solver2->getNumRows();
310
311  int * whichRow = new int[numberRows];
312  // get row copy
313  const CoinPackedMatrix * rowCopy = solver2->getMatrixByRow();
314  const int * column = rowCopy->getIndices();
315  const int * rowLength = rowCopy->getVectorLengths();
316  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
317  const double * rowLower = solver2->getRowLower();
318  const double * rowUpper = solver2->getRowUpper();
319  const double * element = rowCopy->getElements();
320  int iRow,nDelete=0;
321  for (iRow=0;iRow<numberRows;iRow++) {
322    if (rowLower[iRow]<-1.0e20||rowUpper[iRow]>1.0e20) {
323      // take out
324      whichRow[nDelete++]=iRow;
325    }
326  }
327  // leave some rows to avoid empty problem (Gomory does not like)
328  nDelete = CoinMax(CoinMin(nDelete,numberRows-5),0);
329  for (int jRow=0;jRow<nDelete;jRow++) {
330    iRow=whichRow[jRow];
331    int start = rowStart[iRow];
332    stored.addCut(rowLower[iRow],rowUpper[iRow],rowLength[iRow],
333                  column+start,element+start);
334  }
335  /* The following is problem specific.
336     Normally cuts are deleted if slack on cut basic.
337     On some problems you may wish to leave cuts in as long
338     as slack value zero
339  */
340  int numberCuts=stored.sizeRowCuts();
341  for (int iCut=0;iCut<numberCuts;iCut++) {
342    //stored.mutableRowCutPointer(iCut)->setEffectiveness(1.0e50);
343  }
344  solver2->deleteRows(nDelete,whichRow);
345  delete [] whichRow;
346 }
347  CbcModel model(*solver2);
348  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
349  // Set up some cut generators and defaults
350  // Probing first as gets tight bounds on continuous
351
352  CglProbing generator1;
353  generator1.setUsingObjective(true);
354  generator1.setMaxPass(1);
355  generator1.setMaxPassRoot(5);
356  // Number of unsatisfied variables to look at
357  generator1.setMaxProbe(10);
358  generator1.setMaxProbeRoot(1000);
359  // How far to follow the consequences
360  generator1.setMaxLook(50);
361  generator1.setMaxLookRoot(500);
362  // Only look at rows with fewer than this number of elements
363  generator1.setMaxElements(200);
364  generator1.setRowCuts(3);
365
366  CglGomory generator2;
367  // try larger limit
368  generator2.setLimit(300);
369
370  CglKnapsackCover generator3;
371
372  CglRedSplit generator4;
373  // try larger limit
374  generator4.setLimit(200);
375
376  CglClique generator5;
377  generator5.setStarCliqueReport(false);
378  generator5.setRowCliqueReport(false);
379
380  CglMixedIntegerRounding2 mixedGen;
381  CglFlowCover flowGen;
382 
383  // Add in generators
384  // Experiment with -1 and -99 etc
385  // This is just for one particular model
386  model.addCutGenerator(&generator1,-1,"Probing");
387  //model.addCutGenerator(&generator2,-1,"Gomory");
388  model.addCutGenerator(&generator2,1,"Gomory");
389  model.addCutGenerator(&generator3,-1,"Knapsack");
390  // model.addCutGenerator(&generator4,-1,"RedSplit");
391  //model.addCutGenerator(&generator5,-1,"Clique");
392  model.addCutGenerator(&generator5,1,"Clique");
393  model.addCutGenerator(&flowGen,-1,"FlowCover");
394  model.addCutGenerator(&mixedGen,-1,"MixedIntegerRounding");
395  // Add stored cuts (making sure at all depths)
396  model.addCutGenerator(&stored,1,"Stored",true,false,false,-100,1,-1);
397 
398  int numberGenerators = model.numberCutGenerators();
399  int iGenerator;
400  // Say we want timings
401  for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
402    CbcCutGenerator * generator = model.cutGenerator(iGenerator);
403    generator->setTiming(true);
404  }
405  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (model.solver());
406  // go faster stripes
407  if (osiclp) { 
408    if(osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
409      //osiclp->setupForRepeatedUse(2,0);
410      osiclp->setupForRepeatedUse(0,0);
411    }
412    // Don't allow dual stuff
413    osiclp->setSpecialOptions(osiclp->specialOptions()|262144);
414  } 
415  // Uncommenting this should switch off all CBC messages
416  // model.messagesPointer()->setDetailMessages(10,10000,NULL);
417  // No heuristics
418  // Do initial solve to continuous
419  model.initialSolve();
420  /*  You need the next few lines -
421      a) so that cut generator will always be called again if it generated cuts
422      b) it is known that matrix is not enough to define problem so do cuts even
423         if it looks integer feasible at continuous optimum.
424      c) a solution found by strong branching will be ignored.
425      d) don't recompute a solution once found
426  */
427  // Make sure cut generator called correctly (a)
428  iGenerator=numberGenerators-1;
429  model.cutGenerator(iGenerator)->setMustCallAgain(true);
430  // Say cuts needed at continuous (b)
431  OsiBabSolver oddCuts;
432  oddCuts.setSolverType(4);
433  // owing to bug must set after initialSolve
434  model.passInSolverCharacteristics(&oddCuts);
435  // Say no to all solutions by strong branching (c)
436  CbcFeasibilityNoStrong noStrong;
437  model.setProblemFeasibility(noStrong);
438  // Say don't recompute solution d)
439  model.setSpecialOptions(4);
440
441  // Could tune more
442  double objValue = model.solver()->getObjSense()*model.solver()->getObjValue();
443  double minimumDropA=CoinMin(1.0,fabs(objValue)*1.0e-3+1.0e-4);
444  double minimumDrop= fabs(objValue)*1.0e-4+1.0e-4;
445  printf("min drop %g (A %g)\n",minimumDrop,minimumDropA);
446  model.setMinimumDrop(minimumDrop);
447
448  if (model.getNumCols()<500)
449    model.setMaximumCutPassesAtRoot(-100); // always do 100 if possible
450  else if (model.getNumCols()<5000)
451    model.setMaximumCutPassesAtRoot(100); // use minimum drop
452  else
453    model.setMaximumCutPassesAtRoot(20);
454  model.setMaximumCutPasses(10);
455  //model.setMaximumCutPasses(2);
456
457  // Switch off strong branching if wanted
458  // model.setNumberStrong(0);
459  // Do more strong branching if small
460  if (model.getNumCols()<5000)
461    model.setNumberStrong(10);
462  model.setNumberStrong(20);
463  //model.setNumberStrong(5);
464  model.setNumberBeforeTrust(5);
465
466  model.solver()->setIntParam(OsiMaxNumIterationHotStart,100);
467
468  // If time is given then stop after that number of minutes
469  if (minutes>=0.0) {
470    std::cout<<"Stopping after "<<minutes<<" minutes"<<std::endl;
471    model.setDblParam(CbcModel::CbcMaximumSeconds,60.0*minutes);
472  }
473  // Switch off most output
474  if (model.getNumCols()<30000) {
475    model.messageHandler()->setLogLevel(1);
476    //model.solver()->messageHandler()->setLogLevel(0);
477  } else {
478    model.messageHandler()->setLogLevel(2);
479    model.solver()->messageHandler()->setLogLevel(1);
480  }
481  //model.messageHandler()->setLogLevel(2);
482  //model.solver()->messageHandler()->setLogLevel(2);
483  //model.setPrintFrequency(50);
484  //#define DEBUG_CUTS
485#ifdef DEBUG_CUTS
486  // Set up debugger by name (only if no preprocesing)
487  if (!preProcess) {
488    std::string problemName ;
489    model.solver()->getStrParam(OsiProbName,problemName) ;
490    model.solver()->activateRowCutDebugger(problemName.c_str()) ;
491  }
492#endif
493  // Do complete search
494 
495  model.branchAndBound();
496
497  std::cout<<mpsFileName<<" took "<<CoinCpuTime()-time1<<" seconds, "
498           <<model.getNodeCount()<<" nodes with objective "
499           <<model.getObjValue()
500           <<(!model.status() ? " Finished" : " Not finished")
501           <<std::endl;
502
503  // Print more statistics
504  std::cout<<"Cuts at root node changed objective from "<<model.getContinuousObjective()
505           <<" to "<<model.rootObjectiveAfterCuts()<<std::endl;
506
507  for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
508    CbcCutGenerator * generator = model.cutGenerator(iGenerator);
509    std::cout<<generator->cutGeneratorName()<<" was tried "
510             <<generator->numberTimesEntered()<<" times and created "
511             <<generator->numberCutsInTotal()<<" cuts of which "
512             <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
513    if (generator->timing())
514      std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
515    else
516      std::cout<<std::endl;
517  }
518  // Print solution if finished - we can't get names from Osi! - so get from OsiClp
519
520  if (model.getMinimizationObjValue()<1.0e50) {
521    // post process
522    OsiSolverInterface * solver;
523    if (preProcess) {
524      process.postProcess(*model.solver());
525      // Solution now back in solver1
526      solver = & solver1;
527    } else {
528      solver = model.solver();
529    }
530    int numberColumns = solver->getNumCols();
531   
532    const double * solution = solver->getColSolution();
533
534    // Get names from solver1 (as OsiSolverInterface may lose)
535    std::vector<std::string> columnNames = *solver1.getModelPtr()->columnNames();
536   
537    int iColumn;
538    std::cout<<std::setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14);
539   
540    std::cout<<"--------------------------------------"<<std::endl;
541    for (iColumn=0;iColumn<numberColumns;iColumn++) {
542      double value=solution[iColumn];
543      if (fabs(value)>1.0e-7&&solver->isInteger(iColumn)) {
544        std::cout<<std::setw(6)<<iColumn<<" "
545                 <<columnNames[iColumn]<<" "
546                 <<value<<std::endl;
547        solverSave.setColLower(iColumn,value);
548        solverSave.setColUpper(iColumn,value);
549      }
550    }
551    std::cout<<"--------------------------------------"<<std::endl;
552 
553    std::cout<<std::resetiosflags(std::ios::fixed|std::ios::showpoint|std::ios::scientific);
554    solverSave.initialSolve();
555  }
556  return 0;
557}   
Note: See TracBrowser for help on using the repository browser.