source: trunk/Cbc/examples/allCuts.cpp

Last change on this file was 2469, checked in by unxusr, 6 weeks ago

formatting

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