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

Last change on this file since 2496 was 2469, checked in by unxusr, 9 months ago

formatting

  • Property svn:keywords set to Author Date Id Revision
File size: 17.3 KB
RevLine 
[1574]1// $Id: allCuts.cpp 2469 2019-01-06 23:17:46Z forrest $
[953]2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1574]4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
[953]6#include <cassert>
7#include <iomanip>
8
[1898]9#include "CoinPragma.hpp"
[953]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
[1898]39#include "CoinTime.hpp"
[953]40// Need stored cuts
41
42#include "CglStored.hpp"
43
44/** Stored Cut Generator Class */
45class CglStoredUser : public CglStored {
[2469]46
[953]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  */
[2469]58  virtual void generateCuts(const OsiSolverInterface &si, OsiCuts &cs,
59    const CglTreeInfo info = CglTreeInfo()) const;
[953]60  //@}
61
62  /**@name Cut stuff */
63  //@{
[2469]64  OsiRowCut *mutableRowCutPointer(int index)
65  {
66    return cuts_.rowCutPtr(index);
67  }
[953]68  //@}
69
70  /**@name Constructors and destructors */
71  //@{
[2469]72  /// Default constructor
73  CglStoredUser();
[953]74
[2469]75  /// Copy constructor
76  CglStoredUser(const CglStoredUser &rhs);
77
[953]78  /// Clone
[2469]79  virtual CglCutGenerator *clone() const;
[953]80
[2469]81  /// Assignment operator
[953]82  CglStoredUser &
[2469]83  operator=(const CglStoredUser &rhs);
84
85  /// Destructor
86  virtual ~CglStoredUser();
[953]87  //@}
[2469]88
[953]89protected:
[2469]90  // Protected member methods
[953]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
[2469]104//-------------------------------------------------------------------
105void CglStoredUser::generateCuts(const OsiSolverInterface &si, OsiCuts &cs,
106  const CglTreeInfo info) const
[953]107{
108  // Get basic problem information
[2469]109  const double *solution = si.getColSolution();
110  if (info.inTree && info.pass > numberPasses_) {
[953]111    // only continue if integer feasible
[2469]112    int numberColumns = si.getNumCols();
[953]113    int i;
[2469]114    const double *colUpper = si.getColUpper();
115    const double *colLower = si.getColLower();
116    int numberAway = 0;
117    for (i = 0; i < numberColumns; i++) {
[953]118      double value = solution[i];
119      // In case slightly away from bounds
[2469]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++;
[953]124    }
125    if (numberAway)
126      return; // let code branch
127  }
128  int numberRowCuts = cuts_.sizeRowCuts();
[2469]129  for (int i = 0; i < numberRowCuts; i++) {
130    const OsiRowCut *rowCutPointer = cuts_.rowCutPtr(i);
[953]131    double violation = rowCutPointer->violated(solution);
[2469]132    if (violation >= requiredViolation_)
[953]133      cs.insert(*rowCutPointer);
134  }
135}
136
137//-------------------------------------------------------------------
[2469]138// Default Constructor
[953]139//-------------------------------------------------------------------
[2469]140CglStoredUser::CglStoredUser()
141  : CglStored()
142  , numberPasses_(5)
[953]143{
144}
145
146//-------------------------------------------------------------------
[2469]147// Copy constructor
[953]148//-------------------------------------------------------------------
[2469]149CglStoredUser::CglStoredUser(const CglStoredUser &source)
150  : CglStored(source)
151  , numberPasses_(source.numberPasses_)
152{
[953]153}
154
155//-------------------------------------------------------------------
156// Clone
157//-------------------------------------------------------------------
158CglCutGenerator *
159CglStoredUser::clone() const
160{
161  return new CglStoredUser(*this);
162}
163
164//-------------------------------------------------------------------
[2469]165// Destructor
[953]166//-------------------------------------------------------------------
[2469]167CglStoredUser::~CglStoredUser()
[953]168{
169}
170
171//----------------------------------------------------------------
[2469]172// Assignment operator
[953]173//-------------------------------------------------------------------
174CglStoredUser &
[2469]175CglStoredUser::operator=(const CglStoredUser &rhs)
[953]176{
177  if (this != &rhs) {
178    CglStored::operator=(rhs);
[2469]179    numberPasses_ = rhs.numberPasses_;
[953]180  }
181  return *this;
182}
183// Class to disallow strong branching solutions
184#include "CbcFeasibilityBase.hpp"
[2469]185class CbcFeasibilityNoStrong : public CbcFeasibilityBase {
[953]186public:
[2469]187  // Default Constructor
188  CbcFeasibilityNoStrong() {};
[953]189
190  virtual ~CbcFeasibilityNoStrong() {};
[2469]191  // Copy constructor
192  CbcFeasibilityNoStrong(const CbcFeasibilityNoStrong &rhs) {};
[953]193
[2469]194  // Assignment operator
195  CbcFeasibilityNoStrong &operator=(const CbcFeasibilityNoStrong &rhs)
196  {
197    return *this;
198  };
199
[953]200  /// Clone
[2469]201  virtual CbcFeasibilityBase *clone() const
202  {
203    return new CbcFeasibilityNoStrong();
204  };
[953]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  */
[2469]215  virtual int feasible(CbcModel *model, int mode)
216  {
217    return mode;
218  };
[953]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
[2469]231int main(int argc, const char *argv[])
[953]232{
233
234  // Define your favorite OsiSolver
[2469]235
[953]236  OsiClpSolverInterface solver1;
237
238  // Read in model using argv[1]
239  // and assert that it is a clean model
[1464]240  std::string mpsFileName;
[1468]241#if defined(SAMPLEDIR)
[1464]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
[2469]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;
[1898]255  }
[953]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  */
[2469]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;
[953]270      nGoodParam++;
[2469]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;
[953]276          iParam++; // skip time
277        }
278      }
279    }
280  }
[2469]281  if (nGoodParam == 0 && argc == 3 && isdigit(argv[2][0])) {
[953]282    // If time is given then stop after that number of minutes
283    minutes = atof(argv[2]);
[2469]284    if (minutes >= 0.0)
285      nGoodParam = 1;
[953]286  }
[2469]287  if (nGoodParam != argc - 2 && argc >= 2) {
[953]288    printf("Usage <file> [preprocess] [time <minutes>] or <file> <minutes>\n");
289    exit(1);
290  }
291  // Reduce printout
[2469]292  solver1.setHintParam(OsiDoReducePrint, true, OsiHintTry);
[953]293  // See if we want preprocessing
[2469]294  OsiSolverInterface *solver2 = &solver1;
[953]295  CglPreProcess process;
296  // Never do preprocessing until dual tests out as can fix incorrectly
[2469]297  preProcess = false;
[953]298  if (preProcess) {
299    /* Do not try and produce equality cliques and
300       do up to 5 passes */
[2469]301    solver2 = process.preProcess(solver1, false, 5);
[953]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;
[2469]310  {
311    int numberRows = solver2->getNumRows();
[953]312
[2469]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      }
[953]328    }
[2469]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.
[953]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  */
[2469]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;
[953]348  }
349  CbcModel model(*solver2);
[2469]350  model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
[953]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;
[2469]384
[953]385  // Add in generators
386  // Experiment with -1 and -99 etc
387  // This is just for one particular model
[2469]388  model.addCutGenerator(&generator1, -1, "Probing");
[953]389  //model.addCutGenerator(&generator2,-1,"Gomory");
[2469]390  model.addCutGenerator(&generator2, 1, "Gomory");
391  model.addCutGenerator(&generator3, -1, "Knapsack");
[953]392  // model.addCutGenerator(&generator4,-1,"RedSplit");
393  //model.addCutGenerator(&generator5,-1,"Clique");
[2469]394  model.addCutGenerator(&generator5, 1, "Clique");
395  model.addCutGenerator(&flowGen, -1, "FlowCover");
396  model.addCutGenerator(&mixedGen, -1, "MixedIntegerRounding");
[953]397  // Add stored cuts (making sure at all depths)
[2469]398  model.addCutGenerator(&stored, 1, "Stored", true, false, false, -100, 1, -1);
399
[953]400  int numberGenerators = model.numberCutGenerators();
401  int iGenerator;
402  // Say we want timings
[2469]403  for (iGenerator = 0; iGenerator < numberGenerators; iGenerator++) {
404    CbcCutGenerator *generator = model.cutGenerator(iGenerator);
[953]405    generator->setTiming(true);
406  }
[2469]407  OsiClpSolverInterface *osiclp = dynamic_cast< OsiClpSolverInterface * >(model.solver());
[953]408  // go faster stripes
[2469]409  if (osiclp) {
410    if (osiclp->getNumRows() < 300 && osiclp->getNumCols() < 500) {
[953]411      //osiclp->setupForRepeatedUse(2,0);
[2469]412      osiclp->setupForRepeatedUse(0, 0);
[953]413    }
414    // Don't allow dual stuff
[2469]415    osiclp->setSpecialOptions(osiclp->specialOptions() | 262144);
416  }
[953]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)
[2469]430  iGenerator = numberGenerators - 1;
[953]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
[2469]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);
[953]448  model.setMinimumDrop(minimumDrop);
449
[2469]450  if (model.getNumCols() < 500)
[953]451    model.setMaximumCutPassesAtRoot(-100); // always do 100 if possible
[2469]452  else if (model.getNumCols() < 5000)
[953]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
[2469]462  if (model.getNumCols() < 5000)
[953]463    model.setNumberStrong(10);
464  model.setNumberStrong(20);
465  //model.setNumberStrong(5);
466  model.setNumberBeforeTrust(5);
467
[2469]468  model.solver()->setIntParam(OsiMaxNumIterationHotStart, 100);
[953]469
470  // If time is given then stop after that number of minutes
[2469]471  if (minutes >= 0.0) {
472    std::cout << "Stopping after " << minutes << " minutes" << std::endl;
473    model.setDblParam(CbcModel::CbcMaximumSeconds, 60.0 * minutes);
[953]474  }
475  // Switch off most output
[2469]476  if (model.getNumCols() < 30000) {
[953]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) {
[2469]490    std::string problemName;
491    model.solver()->getStrParam(OsiProbName, problemName);
492    model.solver()->activateRowCutDebugger(problemName.c_str());
[953]493  }
494#endif
495  // Do complete search
[2469]496
[953]497  model.branchAndBound();
498
[2469]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;
[953]504
505  // Print more statistics
[2469]506  std::cout << "Cuts at root node changed objective from " << model.getContinuousObjective()
507            << " to " << model.rootObjectiveAfterCuts() << std::endl;
[953]508
[2469]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";
[953]515    if (generator->timing())
[2469]516      std::cout << " ( " << generator->timeInCutGenerator() << " seconds)" << std::endl;
[953]517    else
[2469]518      std::cout << std::endl;
[953]519  }
520  // Print solution if finished - we can't get names from Osi! - so get from OsiClp
521
[2469]522  if (model.getMinimizationObjValue() < 1.0e50) {
[953]523    // post process
[2469]524    OsiSolverInterface *solver;
[953]525    if (preProcess) {
526      process.postProcess(*model.solver());
527      // Solution now back in solver1
[2469]528      solver = &solver1;
[953]529    } else {
530      solver = model.solver();
531    }
532    int numberColumns = solver->getNumCols();
533
[2469]534    const double *solution = solver->getColSolution();
535
[953]536    // Get names from solver1 (as OsiSolverInterface may lose)
[2469]537    std::vector< std::string > columnNames = *solver1.getModelPtr()->columnNames();
538
[953]539    int iColumn;
[2469]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);
[953]551      }
552    }
[2469]553    std::cout << "--------------------------------------" << std::endl;
554
555    std::cout << std::resetiosflags(std::ios::fixed | std::ios::showpoint | std::ios::scientific);
[953]556    solverSave.initialSolve();
557  }
558  return 0;
[2469]559}
Note: See TracBrowser for help on using the repository browser.