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

Last change on this file since 1854 was 1854, checked in by stefan, 6 years ago

fix svn keywords property

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