source: stable/2.1/Cbc/src/unitTestClp.cpp @ 960

Last change on this file since 960 was 960, checked in by forrest, 12 years ago

maximization reporting bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.9 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <cstdio>
5#include <string>
6#include <iostream>
7
8#include "CoinTime.hpp"
9#include "CbcModel.hpp"
10#include "CbcCutGenerator.hpp"
11#include "OsiClpSolverInterface.hpp"
12#include "OsiRowCutDebugger.hpp"
13//#############################################################################
14
15#ifdef NDEBUG
16#undef NDEBUG
17#endif
18
19//#############################################################################
20
21// Display message on stdout and stderr
22void testingMessage( const char * const msg )
23{
24  std::cout <<msg;
25  //cout <<endl <<"*****************************************"
26  //     <<endl <<msg <<endl;
27}
28
29//#############################################################################
30
31static inline bool CbcTestFile(const std::string name)
32{
33  FILE *fp = fopen(name.c_str(),"r");
34  if (fp) {
35    fclose(fp);
36    return true;
37  }
38  return false;
39}
40
41//#############################################################################
42
43bool CbcTestMpsFile(std::string& fname)
44{
45  if (CbcTestFile(fname)) {
46    return true;
47  }
48  if (CbcTestFile(fname+".mps")) {
49    fname += ".mps";
50    return true;
51  }
52  if (CbcTestFile(fname+".MPS")) {
53    fname += ".MPS";
54    return true;
55  }
56#ifdef COIN_HAS_ZLIB
57  if (CbcTestFile(fname+".gz")) {
58    return true;
59  }
60  if (CbcTestFile(fname+".mps.gz")) {
61    fname += ".mps";
62    return true;
63  }
64  if (CbcTestFile(fname+".MPS.gz")) {
65    fname += ".MPS";
66    return true;
67  }
68  if (CbcTestFile(fname+".MPS.GZ")) {
69    fname += ".MPS";
70    return true;
71  }
72#endif
73  return false;
74}
75
76//#############################################################################
77
78int CbcClpUnitTest (const CbcModel & saveModel, std::string& dirMiplib,
79                     bool unitTestOnly)
80{
81  unsigned int m ;
82
83  // Set directory containing miplib data files.
84  std::string test1 = dirMiplib +"p0033";
85  // See if files exist
86  bool doTest=CbcTestMpsFile(test1);
87
88  if (!doTest) {
89    printf("Not doing miplib run as can't find mps files - ? .gz without libz\n");
90    return -1;
91  }
92  /*
93    Vectors to hold test problem names and characteristics. The objective value
94    after optimization (objValue) must agree to the specified tolerance
95    (objValueTol).
96  */
97  std::vector<std::string> mpsName ;
98  std::vector<int> nRows ;
99  std::vector<int> nCols ;
100  std::vector<double> objValueC ;
101  std::vector<double> objValue ;
102  std::vector<int> strategy ;
103  /*
104    And a macro to make the vector creation marginally readable.
105  */
106#define PUSH_MPS(zz_mpsName_zz,\
107                 zz_nRows_zz,zz_nCols_zz,zz_objValue_zz,zz_objValueC_zz, \
108                 zz_strategy_zz) \
109  mpsName.push_back(zz_mpsName_zz) ; \
110  nRows.push_back(zz_nRows_zz) ; \
111  nCols.push_back(zz_nCols_zz) ; \
112  objValueC.push_back(zz_objValueC_zz) ; \
113  strategy.push_back(zz_strategy_zz) ; \
114  objValue.push_back(zz_objValue_zz) ;
115
116  if (unitTestOnly) {
117    PUSH_MPS("p0033",16,33,3089,2520.57,7);
118    PUSH_MPS("p0201",133,201,7615,6875.0,7);
119  } else {
120    /*
121      Load up the problem vector. Note that the row counts here include the
122      objective function.
123    */
124    // 0 for no test, 1 for some, 2 for many, 3 for all
125    //PUSH_MPS("blend2",274,353,7.598985,6.9156751140,7);
126    //PUSH_MPS("egout",98,141,568.101,149.589,7);
127#define HOWMANY 2
128#if HOWMANY
129#if HOWMANY>1
130    PUSH_MPS("10teams",230,2025,924,917,7);
131#endif
132    PUSH_MPS("air03",124,10757,340160,338864.25,7);
133#if HOWMANY>2
134    PUSH_MPS("air04",823,8904,56137,55535.436,8);
135    PUSH_MPS("air05",426,7195,26374,25877.609,8);
136#endif
137  //    PUSH_MPS("arki001",1048,1388,7580813.0459,7579599.80787,7);
138    PUSH_MPS("bell3a",123,133,878430.32,862578.64,7);
139#if HOWMANY>1
140    PUSH_MPS("bell5",91,104,8966406.49,8608417.95,7);
141#endif
142    PUSH_MPS("blend2",274,353,7.598985,6.9156751140,7);
143#if HOWMANY>1
144    PUSH_MPS("cap6000",2176,6000,-2451377,-2451537.325,7);
145#endif
146    //    PUSH_MPS("dano3mip",3202,13873,728.1111,576.23162474,7);
147    //PUSH_MPS("danoint",664,521,65.67,62.637280418,7);
148    PUSH_MPS("dcmulti",290,548,188182,183975.5397,7);
149    PUSH_MPS("dsbmip",1182,1886,-305.19817501,-305.19817501,7);
150    PUSH_MPS("egout",98,141,568.101,149.589,7);
151    PUSH_MPS("enigma",21,100,0.0,0.0,7);
152#if HOWMANY>3
153    PUSH_MPS("fast0507",507,63009,174,172.14556668,7);
154#endif
155    PUSH_MPS("fiber",363,1298,405935.18000,156082.51759,7);
156#if HOWMANY>1
157    PUSH_MPS("fixnet6",478,878,3983,1200.88,7);
158#endif
159    PUSH_MPS("flugpl",18,18,1201500,1167185.7,7);
160    PUSH_MPS("gen",780,870,112313,112130.0,7);
161#if HOWMANY>1
162    PUSH_MPS("gesa2",1392,1224,25779856.372,25476489.678,7);
163    PUSH_MPS("gesa2_o",1248,1224,25779856.372,25476489.678,7);
164#endif
165    PUSH_MPS("gesa3",1368,1152,27991042.648,27833632.451,7);
166    PUSH_MPS("gesa3_o",1224,1152,27991042.648,27833632.451,7);
167    PUSH_MPS("gt2",29,188,21166.000,13460.233074,7);
168#if HOWMANY>2
169    PUSH_MPS("harp2",112,2993,-73899798.00,-74353341.502,7);
170#endif
171    PUSH_MPS("khb05250",101,1350,106940226,95919464.0,7);
172#if HOWMANY>1
173    PUSH_MPS("l152lav",97,1989,4722,4656.36,7);
174#endif
175    PUSH_MPS("lseu",28,89,1120,834.68,7);
176    PUSH_MPS("misc03",96,160,3360,1910.,7);
177    PUSH_MPS("misc06",820,1808,12850.8607,12841.6,7);
178#if HOWMANY>1
179    PUSH_MPS("misc07",212,260,2810,1415.0,7);
180    PUSH_MPS("mitre",2054,10724,115155,114740.5184,7);
181#endif
182    PUSH_MPS("mod008",6,319,307,290.9,7);
183    PUSH_MPS("mod010",146,2655,6548,6532.08,7);
184#if HOWMANY>2
185    PUSH_MPS("mod011",4480,10958,-54558535,-62121982.55,7);
186    PUSH_MPS("modglob",291,422,20740508,20430947.,7);
187#endif
188#if HOWMANY>3
189    PUSH_MPS("noswot",182,128,-43,-43.0,7);
190#endif
191#if HOWMANY>1
192    PUSH_MPS("nw04",36,87482,16862,16310.66667,7);
193#endif
194    PUSH_MPS("p0033",16,33,3089,2520.57,7);
195    PUSH_MPS("p0201",133,201,7615,6875.0,7);
196    PUSH_MPS("p0282",241,282,258411,176867.50,7);
197    PUSH_MPS("p0548",176,548,8691,315.29,7);
198    PUSH_MPS("p2756",755,2756,3124,2688.75,7);
199#if HOWMANY>2
200    PUSH_MPS("pk1",45,86,11.0,0.0,7);
201#endif
202#if HOWMANY>1
203    PUSH_MPS("pp08a",136,240,7350.0,2748.3452381,7);
204    PUSH_MPS("pp08aCUTS",246,240,7350.0,5480.6061563,7);
205#endif
206#if HOWMANY>2
207    PUSH_MPS("qiu",1192,840,-132.873137,-931.638857,7);
208#endif
209    PUSH_MPS("qnet1",503,1541,16029.692681,14274.102667,7);
210    PUSH_MPS("qnet1_o",456,1541,16029.692681,12095.571667,7);
211    PUSH_MPS("rentacar",6803,9557,30356761,28806137.644,7);
212    PUSH_MPS("rgn",24,180,82.1999,48.7999,7);
213#if HOWMANY>3
214    PUSH_MPS("rout",291,556,1077.56,981.86428571,7);
215    PUSH_MPS("set1ch",492,712,54537.75,32007.73,7);
216#endif
217    //    PUSH_MPS("seymour",4944,1372,423,403.84647413,7);
218    PUSH_MPS("stein27",118,27,18,13.0,7);
219#if HOWMANY>1
220    PUSH_MPS("stein45",331,45,30,22.0,7);
221#endif
222    PUSH_MPS("vpm1",234,378,20,15.4167,7);
223    PUSH_MPS("vpm2",234,378,13.75,9.8892645972,7);
224#endif
225  }
226#undef PUSH_MPS
227   
228  int numProbSolved = 0;
229  double timeTaken=0.0;
230  int numberFailures=0;
231 
232  /*
233    Open the main loop to step through the MPS problems.
234  */
235  for (m = 0 ; m < mpsName.size() ; m++) {
236    std::cout << "  processing mps file: " << mpsName[m] 
237              << " (" << m+1 << " out of " << mpsName.size() << ")\n";
238    /*
239      Stage 1: Read the MPS
240      and make sure the size of the constraint matrix is correct.
241    */
242    CbcModel * model = new CbcModel(saveModel);
243     
244    std::string fn = dirMiplib+mpsName[m] ;
245    if (!CbcTestMpsFile(fn)) {
246      std::cout << "ERROR: Cannot find MPS file " << fn << "\n";
247      continue;
248    }
249    model->solver()->readMps(fn.c_str(),"") ;
250    assert(model->getNumRows() == nRows[m]) ;
251    assert(model->getNumCols() == nCols[m]) ;
252
253    /*
254      Stage 2: Call solver to solve the problem.
255      then check the return code and objective.
256    */
257
258    double startTime = CoinCpuTime()+CoinCpuTimeJustChildren();
259    if (model->getMaximumNodes()>200000) {
260      model->setMaximumNodes(200000);
261    }
262    OsiClpSolverInterface * si =
263      dynamic_cast<OsiClpSolverInterface *>(model->solver()) ;
264    assert (si != NULL);
265    // get clp itself
266    ClpSimplex * modelC = si->getModelPtr();
267    modelC->tightenPrimalBounds(0.0,0,true);
268    model->initialSolve();
269    if (modelC->dualBound()==1.0e10) {
270      // user did not set - so modify
271      // get largest scaled away from bound
272      double largest=1.0e-12;
273      int numberRows = modelC->numberRows();
274      const double * rowPrimal = modelC->primalRowSolution();
275      const double * rowLower = modelC->rowLower();
276      const double * rowUpper = modelC->rowUpper();
277      const double * rowScale = modelC->rowScale();
278      int iRow;
279      for (iRow=0;iRow<numberRows;iRow++) {
280        double value = rowPrimal[iRow];
281        double above = value-rowLower[iRow];
282        double below = rowUpper[iRow]-value;
283        if (rowScale) {
284          double multiplier = rowScale[iRow];
285          above *= multiplier;
286          below *= multiplier;
287        }
288        if (above<1.0e12) {
289          largest = CoinMax(largest,above);
290        }
291        if (below<1.0e12) {
292          largest = CoinMax(largest,below);
293        }
294      }
295     
296      int numberColumns = modelC->numberColumns();
297      const double * columnPrimal = modelC->primalColumnSolution();
298      const double * columnLower = modelC->columnLower();
299      const double * columnUpper = modelC->columnUpper();
300      const double * columnScale = modelC->columnScale();
301      int iColumn;
302      for (iColumn=0;iColumn<numberColumns;iColumn++) {
303        double value = columnPrimal[iColumn];
304        double above = value-columnLower[iColumn];
305        double below = columnUpper[iColumn]-value;
306        if (columnScale) {
307          double multiplier = 1.0/columnScale[iColumn];
308          above *= multiplier;
309          below *= multiplier;
310        }
311        if (above<1.0e12) {
312          largest = CoinMax(largest,above);
313        }
314        if (below<1.0e12) {
315          largest = CoinMax(largest,below);
316        }
317      }
318      //std::cout<<"Largest (scaled) away from bound "<<largest<<std::endl;
319      modelC->setDualBound(CoinMax(1.0001e8,
320                                   CoinMin(1000.0*largest,1.00001e10)));
321    }
322    model->setMinimumDrop(min(5.0e-2,
323                              fabs(model->getMinimizationObjValue())*1.0e-3+1.0e-4));
324    if (model->getNumCols()<500) {
325      model->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
326    } else if (model->getNumCols()<5000) {
327      model->setMaximumCutPassesAtRoot(100); // use minimum drop
328    } else {
329      model->setMaximumCutPassesAtRoot(20);
330    }
331    // If defaults then increase trust for small models
332    if (model->numberStrong()==5&&model->numberBeforeTrust()==10) {
333      int numberColumns = model->getNumCols();
334      if (numberColumns<=50) {
335        model->setNumberBeforeTrust(1000);
336      } else if (numberColumns<=100) {
337        model->setNumberBeforeTrust(100);
338      } else if (numberColumns<=300) {
339        model->setNumberBeforeTrust(50);
340      }
341    }
342    if (model->getNumCols()==-353) {
343      // blend2
344      std::string problemName ;
345      model->solver()->getStrParam(OsiProbName,problemName) ;
346      model->solver()->activateRowCutDebugger(problemName.c_str()) ;
347    }
348    setCutAndHeuristicOptions(*model);
349    model->branchAndBound();
350     
351    double timeOfSolution = CoinCpuTime()+CoinCpuTimeJustChildren()-startTime;
352    // Print more statistics
353    std::cout<<"Cuts at root node changed objective from "<<model->getContinuousObjective()
354             <<" to "<<model->rootObjectiveAfterCuts()<<std::endl;
355    int numberGenerators = model->numberCutGenerators();
356    for (int iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
357      CbcCutGenerator * generator = model->cutGenerator(iGenerator);
358      std::cout<<generator->cutGeneratorName()<<" was tried "
359               <<generator->numberTimesEntered()<<" times and created "
360               <<generator->numberCutsInTotal()<<" cuts of which "
361               <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
362      if (generator->timing())
363        std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
364      else
365        std::cout<<std::endl;
366    }
367    if (!model->status()) { 
368      double soln = model->getObjValue();       
369      CoinRelFltEq eq(1.0e-3) ;
370      if (eq(soln,objValue[m])) { 
371        std::cout
372          <<"cbc_clp"<<" "
373          << soln << " = " << objValue[m] << " ; okay";
374        numProbSolved++;
375      } else  { 
376        std::cout <<"cbc_clp" <<" " <<soln << " != " <<objValue[m]
377                  << "; error=" << fabs(objValue[m] - soln); 
378        numberFailures++;
379        //#ifdef COIN_DEVELOP
380        //abort();
381        //#endif
382      }
383    } else {
384      std::cout << "error; too many nodes" ;
385    }
386    std::cout<<" - took " <<timeOfSolution<<" seconds.("<<
387      model->getNodeCount()<<" / "<<model->getIterationCount()<<
388      " )"<<std::endl;
389    timeTaken += timeOfSolution;
390    delete model;
391  }
392  int returnCode=0;
393  std::cout
394    <<"cbc_clp" 
395    <<" solved " 
396    <<numProbSolved
397    <<" out of "
398    <<objValue.size();
399  int numberOnNodes = objValue.size()-numProbSolved-numberFailures;
400  if (numberFailures||numberOnNodes) {
401    if (numberOnNodes) {
402      std::cout<<" ("<<numberOnNodes<<" stopped on nodes)";
403      returnCode = numberOnNodes;
404    }
405    if (numberFailures) {
406      std::cout<<" ("<<numberFailures<<" gave bad answer!)";
407      returnCode += 100*numberFailures;
408    }
409  }
410  std::cout<<" and took "
411    <<timeTaken
412    <<" seconds."
413    <<std::endl;
414  return returnCode;
415}
Note: See TracBrowser for help on using the repository browser.