← Back to branch summary

~ubuntu-branches/ubuntu/natty/coinor-cbc/natty

« back to all changes in this revision

Viewing changes to Cbc/examples/repeat.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Soeren Sonnenburg
  • Date: 2009-08-26 16:57:29 UTC
  • Revision ID: james.westby@ubuntu.com-20090826165729-ybrezxdybg0hd1u6
Tags: upstream-2.3.0
Import upstream version 2.3.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Cbc/examples/repeat.cpp
 
1
// Copyright (C) 2005, 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 "CbcStrategy.hpp"
 
16
#include "CbcBranchUser.hpp"
 
17
#include "CbcCompareUser.hpp"
 
18
#include "CbcCutGenerator.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
// Heuristics
 
35
 
 
36
#include "CbcHeuristic.hpp"
 
37
 
 
38
#include  "CoinTime.hpp"
 
39
 
 
40
//#############################################################################
 
41
 
 
42
 
 
43
/************************************************************************
 
44
 
 
45
This main program reads in an integer model from an mps file.
 
46
 
 
47
It then sets up some Cgl cut generators and calls branch and cut.
 
48
 
 
49
This is designed to show two things :-
 
50
 
 
51
1) Use of CbcStrategy to do preprocessing - this should be cleaner than old way
 
52
2) Looping round modifying solver and redoing branch and cut.  This uses resetToReferenceCopy
 
53
   which resets stuff like the cutoff and any memory of previous branch and cut.
 
54
 
 
55
************************************************************************/
 
56
 
 
57
 
 
58
int main (int argc, const char *argv[])
 
59
{
 
60
 
 
61
  // Define your favorite OsiSolver
 
62
  
 
63
  OsiClpSolverInterface solver1;
 
64
 
 
65
  // Read in model using argv[1]
 
66
  // and assert that it is a clean model
 
67
  std::string mpsFileName = "../../Data/Sample/p0033.mps";
 
68
  if (argc>=2) mpsFileName = argv[1];
 
69
  int numMpsReadErrors = solver1.readMps(mpsFileName.c_str(),"");
 
70
  assert(numMpsReadErrors==0);
 
71
  double time1 = CoinCpuTime();
 
72
 
 
73
  /* Options are:
 
74
     preprocess to do preprocessing
 
75
     time in minutes
 
76
     if 2 parameters and numeric taken as time
 
77
  */
 
78
  bool preProcess=false;
 
79
  double minutes=-1.0;
 
80
  int nGoodParam=0;
 
81
  for (int iParam=2; iParam<argc;iParam++) {
 
82
    if (!strcmp(argv[iParam],"preprocess")) {
 
83
      preProcess=true;
 
84
      nGoodParam++;
 
85
    } else if (!strcmp(argv[iParam],"time")) {
 
86
      if (iParam+1<argc&&isdigit(argv[iParam+1][0])) {
 
87
        minutes=atof(argv[iParam+1]);
 
88
        if (minutes>=0.0) {
 
89
          nGoodParam+=2;
 
90
          iParam++; // skip time
 
91
        }
 
92
      }
 
93
    }
 
94
  }
 
95
  if (nGoodParam==0&&argc==3&&isdigit(argv[2][0])) {
 
96
    // If time is given then stop after that number of minutes
 
97
    minutes = atof(argv[2]);
 
98
    if (minutes>=0.0) 
 
99
      nGoodParam=1;
 
100
  }
 
101
  if (nGoodParam!=argc-2&&argc>=2) {
 
102
    printf("Usage <file> [preprocess] [time <minutes>] or <file> <minutes>\n");
 
103
    exit(1);
 
104
  }
 
105
  //solver1.getModelPtr()->setLogLevel(0);
 
106
  solver1.messageHandler()->setLogLevel(0);
 
107
  solver1.initialSolve();
 
108
  // Reduce printout
 
109
  solver1.setHintParam(OsiDoReducePrint,true,OsiHintTry);
 
110
  CbcModel model(solver1);
 
111
  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
 
112
  // Set up some cut generators and defaults
 
113
  // Probing first as gets tight bounds on continuous
 
114
 
 
115
  CglProbing generator1;
 
116
  generator1.setUsingObjective(true);
 
117
  generator1.setMaxPass(1);
 
118
  generator1.setMaxPassRoot(5);
 
119
  // Number of unsatisfied variables to look at
 
120
  generator1.setMaxProbe(10);
 
121
  generator1.setMaxProbeRoot(1000);
 
122
  // How far to follow the consequences
 
123
  generator1.setMaxLook(50);
 
124
  generator1.setMaxLookRoot(500);
 
125
  // Only look at rows with fewer than this number of elements
 
126
  generator1.setMaxElements(200);
 
127
  generator1.setRowCuts(3);
 
128
 
 
129
  CglGomory generator2;
 
130
  // try larger limit
 
131
  generator2.setLimit(300);
 
132
 
 
133
  CglKnapsackCover generator3;
 
134
 
 
135
  CglRedSplit generator4;
 
136
  // try larger limit
 
137
  generator4.setLimit(200);
 
138
 
 
139
  CglClique generator5;
 
140
  generator5.setStarCliqueReport(false);
 
141
  generator5.setRowCliqueReport(false);
 
142
 
 
143
  CglMixedIntegerRounding2 mixedGen;
 
144
  CglFlowCover flowGen;
 
145
  
 
146
  // Add in generators
 
147
  // Experiment with -1 and -99 etc
 
148
  model.addCutGenerator(&generator1,-1,"Probing");
 
149
  model.addCutGenerator(&generator2,-1,"Gomory");
 
150
  model.addCutGenerator(&generator3,-1,"Knapsack");
 
151
  // model.addCutGenerator(&generator4,-1,"RedSplit");
 
152
  model.addCutGenerator(&generator5,-1,"Clique");
 
153
  model.addCutGenerator(&flowGen,-1,"FlowCover");
 
154
  model.addCutGenerator(&mixedGen,-1,"MixedIntegerRounding");
 
155
  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (model.solver());
 
156
  // go faster stripes
 
157
  if (osiclp) {
 
158
    // Turn this off if you get problems
 
159
    // Used to be automatically set
 
160
    osiclp->setSpecialOptions(128);
 
161
    if(osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
 
162
      //osiclp->setupForRepeatedUse(2,1);
 
163
      osiclp->setupForRepeatedUse(0,1);
 
164
    }
 
165
  } 
 
166
  // Uncommenting this should switch off most CBC messages
 
167
  //model.messagesPointer()->setDetailMessages(10,5,5000);
 
168
  // Allow rounding heuristic
 
169
 
 
170
  CbcRounding heuristic1(model);
 
171
  model.addHeuristic(&heuristic1);
 
172
 
 
173
  // And local search when new solution found
 
174
 
 
175
  CbcHeuristicLocal heuristic2(model);
 
176
  model.addHeuristic(&heuristic2);
 
177
 
 
178
  // Redundant definition of default branching (as Default == User)
 
179
  CbcBranchUserDecision branch;
 
180
  model.setBranchingMethod(&branch);
 
181
 
 
182
  // Definition of node choice
 
183
  CbcCompareUser compare;
 
184
  model.setNodeComparison(compare);
 
185
 
 
186
  // Do initial solve to continuous
 
187
  model.initialSolve();
 
188
 
 
189
  // Could tune more
 
190
  double objValue = model.solver()->getObjSense()*model.solver()->getObjValue();
 
191
  double minimumDropA=CoinMin(1.0,fabs(objValue)*1.0e-3+1.0e-4);
 
192
  double minimumDrop= fabs(objValue)*1.0e-4+1.0e-4;
 
193
  printf("min drop %g (A %g)\n",minimumDrop,minimumDropA);
 
194
  model.setMinimumDrop(minimumDrop);
 
195
 
 
196
  if (model.getNumCols()<500)
 
197
    model.setMaximumCutPassesAtRoot(-100); // always do 100 if possible
 
198
  else if (model.getNumCols()<5000)
 
199
    model.setMaximumCutPassesAtRoot(100); // use minimum drop
 
200
  else
 
201
    model.setMaximumCutPassesAtRoot(20);
 
202
  model.setMaximumCutPasses(10);
 
203
  //model.setMaximumCutPasses(2);
 
204
 
 
205
  // Switch off strong branching if wanted
 
206
  // model.setNumberStrong(0);
 
207
  // Do more strong branching if small
 
208
  if (model.getNumCols()<5000)
 
209
    model.setNumberStrong(10);
 
210
  model.setNumberStrong(20);
 
211
  //model.setNumberStrong(5);
 
212
  model.setNumberBeforeTrust(5);
 
213
  //model.setSizeMiniTree(2);
 
214
 
 
215
  model.solver()->setIntParam(OsiMaxNumIterationHotStart,100);
 
216
 
 
217
  // If time is given then stop after that number of minutes
 
218
  if (minutes>=0.0) {
 
219
    std::cout<<"Stopping after "<<minutes<<" minutes"<<std::endl;
 
220
    model.setDblParam(CbcModel::CbcMaximumSeconds,60.0*minutes);
 
221
  }
 
222
  // Switch off most output
 
223
  if (model.getNumCols()<3000) {
 
224
    model.messageHandler()->setLogLevel(1);
 
225
    //model.solver()->messageHandler()->setLogLevel(0);
 
226
  } else {
 
227
    model.messageHandler()->setLogLevel(2);
 
228
    model.solver()->messageHandler()->setLogLevel(1);
 
229
  }
 
230
  // Default strategy will leave cut generators as they exist already
 
231
  // so cutsOnlyAtRoot (1) ignored
 
232
  // numberStrong (2) is 5 (default)
 
233
  // numberBeforeTrust (3) is 5 (default is 0)
 
234
  // printLevel (4) defaults (0)
 
235
  CbcStrategyDefault strategy(true,5,5);
 
236
  // Set up pre-processing to find sos if wanted
 
237
  if (preProcess)
 
238
    strategy.setupPreProcessing(2);
 
239
  model.setStrategy(strategy);
 
240
 
 
241
  // Go round adding cuts to cutoff last solution
 
242
  // Stop after finding 20 best solutions
 
243
  for (int iPass=0;iPass<20;iPass++) {
 
244
    time1 = CoinCpuTime();
 
245
    // Do complete search
 
246
    model.branchAndBound();
 
247
    
 
248
    std::cout<<mpsFileName<<" took "<<CoinCpuTime()-time1<<" seconds, "
 
249
             <<model.getNodeCount()<<" nodes with objective "
 
250
             <<model.getObjValue()
 
251
             <<(!model.status() ? " Finished" : " Not finished")
 
252
             <<std::endl;
 
253
    // Stop if infeasible
 
254
    if (model.isProvenInfeasible())
 
255
      break;
 
256
    // Print solution if finished - we can't get names from Osi! - so get from OsiClp
 
257
    
 
258
    assert (model.getMinimizationObjValue()<1.0e50);
 
259
    OsiSolverInterface * solver = model.solver();
 
260
    int numberColumns = solver->getNumCols();
 
261
    
 
262
    const double * solution = model.bestSolution();
 
263
    //const double * lower = solver->getColLower();
 
264
    //const double * upper = solver->getColUpper();
 
265
 
 
266
    // Get names from solver1 (as OsiSolverInterface may lose)
 
267
    std::vector<std::string> columnNames = *solver1.getModelPtr()->columnNames();
 
268
    
 
269
    int iColumn;
 
270
    std::cout<<std::setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14);
 
271
    
 
272
    std::cout<<"--------------------------------------"<<std::endl;
 
273
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
 
274
      double value=solution[iColumn];
 
275
      if (fabs(value)>1.0e-7&&solver->isInteger(iColumn)) 
 
276
        std::cout<<std::setw(6)<<iColumn<<" "
 
277
                 <<columnNames[iColumn]<<" "
 
278
                 <<value
 
279
          //<<" "<<lower[iColumn]<<" "<<upper[iColumn]
 
280
                 <<std::endl;
 
281
    }
 
282
    std::cout<<"--------------------------------------"<<std::endl;
 
283
  
 
284
    std::cout<<std::resetiosflags(std::ios::fixed|std::ios::showpoint|std::ios::scientific);
 
285
    /* Now add cut to reference copy.
 
286
       resetting to reference copy also gets rid of best solution so we
 
287
       should either save best solution, reset, add cut OR
 
288
       add cut to reference copy then reset - this is doing latter
 
289
    */
 
290
    OsiSolverInterface * refSolver = model.referenceSolver();
 
291
    const double * bestSolution = model.bestSolution();
 
292
    const double * originalLower = refSolver->getColLower();
 
293
    const double * originalUpper = refSolver->getColUpper();
 
294
    CoinPackedVector cut;
 
295
    double rhs = 1.0;
 
296
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
 
297
      double value=bestSolution[iColumn];
 
298
      if (solver->isInteger(iColumn)) {
 
299
        // only works for 0-1 variables
 
300
        assert (originalLower[iColumn]==0.0&&
 
301
                originalUpper[iColumn]==1.0);
 
302
        // double check integer
 
303
        assert (fabs(floor(value+0.5)-value)<1.0e-5);
 
304
        if (value>0.5) {
 
305
          // at 1.0
 
306
          cut.insert(iColumn,-1.0);
 
307
          rhs -= 1.0;
 
308
        } else {
 
309
          // at 0.0
 
310
          cut.insert(iColumn,1.0);
 
311
        }
 
312
      }
 
313
    }
 
314
    // now add cut
 
315
    refSolver->addRow(cut,rhs,COIN_DBL_MAX);
 
316
    model.resetToReferenceSolver();
 
317
  }
 
318
  return 0;
 
319
}