source: trunk/Cbc/examples/allCuts.cpp

Last change on this file was 1898, checked in by stefan, 5 years ago

fixup examples

  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1// $Id: allCuts.cpp 1898 2013-04-09 18:06:04Z forrest $
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#include <cassert>
7#include <iomanip>
8
9
10#include "CoinPragma.hpp"
11// For Branch and bound
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcBranchUser.hpp"
15#include "CbcCompareUser.hpp"
16#include "CbcCutGenerator.hpp"
17#include "CbcStrategy.hpp"
18#include "CbcHeuristicLocal.hpp"
19#include "OsiClpSolverInterface.hpp"
20
21// Cuts
22
23#include "CglGomory.hpp"
24#include "CglProbing.hpp"
25#include "CglKnapsackCover.hpp"
26#include "CglRedSplit.hpp"
27#include "CglClique.hpp"
28#include "CglFlowCover.hpp"
29#include "CglMixedIntegerRounding2.hpp"
30// Preprocessing
31#include "CglPreProcess.hpp"
32
33// For saying about solution validity
34#include "OsiAuxInfo.hpp"
35
36// Heuristics (But any will have to be special)
37
38#include "CbcHeuristic.hpp"
39
40#include "CoinTime.hpp"
41// Need stored cuts
42
43#include "CglStored.hpp"
44
45/** Stored Cut Generator Class */
46class CglStoredUser : public CglStored {
47 
48public:
49   
50 
51  /**@name Generate Cuts */
52  //@{
53  /** Generate Mixed Integer Stored cuts for the model of the
54      solver interface, si.
55
56      Insert the generated cuts into OsiCut, cs.
57
58      This generator just looks at previously stored cuts
59      and inserts any that are violated by enough
60  */
61  virtual void generateCuts( const OsiSolverInterface & si, OsiCuts & cs,
62                             const CglTreeInfo info = CglTreeInfo()) const;
63  //@}
64
65  /**@name Cut stuff */
66  //@{
67  OsiRowCut * mutableRowCutPointer(int index)
68  { return cuts_.rowCutPtr(index);}
69  //@}
70
71  /**@name Constructors and destructors */
72  //@{
73  /// Default constructor
74  CglStoredUser ();
75 
76  /// Copy constructor
77  CglStoredUser (const CglStoredUser & rhs);
78
79  /// Clone
80  virtual CglCutGenerator * clone() const;
81
82  /// Assignment operator
83  CglStoredUser &
84    operator=(const CglStoredUser& rhs);
85 
86  /// Destructor
87  virtual
88    ~CglStoredUser ();
89  //@}
90     
91protected:
92 
93 // Protected member methods
94
95  // Protected member data
96
97  /**@name Protected member data */
98  //@{
99  /** Don't add any more cuts after this number passes (per node)
100      unless looks integer feasible
101  */
102  int numberPasses_;
103  //@}
104};
105//-------------------------------------------------------------------
106// Generate Stored cuts
107//-------------------------------------------------------------------
108void 
109CglStoredUser::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
110                             const CglTreeInfo info) const
111{
112  // Get basic problem information
113  const double * solution = si.getColSolution();
114  if (info.inTree&&info.pass>numberPasses_) {
115    // only continue if integer feasible
116    int numberColumns=si.getNumCols(); 
117    int i;
118    const double * colUpper = si.getColUpper();
119    const double * colLower = si.getColLower();
120    int numberAway=0;
121    for (i=0;i<numberColumns;i++) {
122      double value = solution[i];
123      // In case slightly away from bounds
124      value = CoinMax(colLower[i],value);
125      value = CoinMin(colUpper[i],value);
126      if (si.isInteger(i)&&fabs(value-fabs(value+0.5))>1.0e-5) 
127        numberAway++;
128    }
129    if (numberAway)
130      return; // let code branch
131  }
132  int numberRowCuts = cuts_.sizeRowCuts();
133  for (int i=0;i<numberRowCuts;i++) {
134    const OsiRowCut * rowCutPointer = cuts_.rowCutPtr(i);
135    double violation = rowCutPointer->violated(solution);
136    if (violation>=requiredViolation_)
137      cs.insert(*rowCutPointer);
138  }
139}
140
141//-------------------------------------------------------------------
142// Default Constructor
143//-------------------------------------------------------------------
144CglStoredUser::CglStoredUser ()
145:
146  CglStored(),
147  numberPasses_(5)
148{
149}
150
151//-------------------------------------------------------------------
152// Copy constructor
153//-------------------------------------------------------------------
154CglStoredUser::CglStoredUser (const CglStoredUser & source) :
155  CglStored(source),
156  numberPasses_(source.numberPasses_)
157{ 
158}
159
160//-------------------------------------------------------------------
161// Clone
162//-------------------------------------------------------------------
163CglCutGenerator *
164CglStoredUser::clone() const
165{
166  return new CglStoredUser(*this);
167}
168
169//-------------------------------------------------------------------
170// Destructor
171//-------------------------------------------------------------------
172CglStoredUser::~CglStoredUser ()
173{
174}
175
176//----------------------------------------------------------------
177// Assignment operator
178//-------------------------------------------------------------------
179CglStoredUser &
180CglStoredUser::operator=(const CglStoredUser& rhs)
181{
182  if (this != &rhs) {
183    CglStored::operator=(rhs);
184    numberPasses_=rhs.numberPasses_;
185  }
186  return *this;
187}
188// Class to disallow strong branching solutions
189#include "CbcFeasibilityBase.hpp"
190class CbcFeasibilityNoStrong : public CbcFeasibilityBase{
191public:
192  // Default Constructor
193  CbcFeasibilityNoStrong () {};
194
195  virtual ~CbcFeasibilityNoStrong() {};
196  // Copy constructor
197  CbcFeasibilityNoStrong ( const CbcFeasibilityNoStrong &rhs) {};
198   
199  // Assignment operator
200  CbcFeasibilityNoStrong & operator=( const CbcFeasibilityNoStrong& rhs)
201  { return * this;};
202
203  /// Clone
204  virtual CbcFeasibilityBase * clone() const
205  { return new CbcFeasibilityNoStrong();};
206
207  /**
208     On input mode:
209     0 - called after a solve but before any cuts
210     -1 - called after strong branching
211     Returns :
212     0 - no opinion
213     -1 pretend infeasible
214     1 pretend integer solution
215  */
216  virtual int feasible(CbcModel * model, int mode) 
217  {return mode;};
218};
219
220
221//#############################################################################
222
223
224/************************************************************************
225
226This main program reads in an integer model from an mps file.
227
228It makes L or G rows into cuts
229
230************************************************************************/
231
232int main (int argc, const char *argv[])
233{
234
235  // Define your favorite OsiSolver
236 
237  OsiClpSolverInterface solver1;
238
239  // Read in model using argv[1]
240  // and assert that it is a clean model
241  std::string mpsFileName;
242#if defined(SAMPLEDIR)
243  mpsFileName = SAMPLEDIR "/p0033.mps";
244#else
245  if (argc < 2) {
246    fprintf(stderr, "Do not know where to find sample MPS files.\n");
247    exit(1);
248  }
249#endif
250  if (argc>=2) mpsFileName = argv[1];
251  int numMpsReadErrors = solver1.readMps(mpsFileName.c_str(),"");
252  if( numMpsReadErrors != 0 )
253  {
254     printf("%d errors reading MPS file\n", numMpsReadErrors);
255     return numMpsReadErrors;
256  }
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.