source: trunk/Cbc/src/unitTestClp.cpp @ 935

Last change on this file since 935 was 935, checked in by forrest, 13 years ago

dd message to say unit test passed

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