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

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

fix bug in unitTest

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