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

Last change on this file since 2370 was 2370, checked in by forrest, 3 years ago

changes for rays

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.3 KB
Line 
1/* $Id: unitTestClp.cpp 2370 2018-03-31 16:45:28Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7#include <cstdio>
8#include <string>
9#include <iostream>
10#include <iomanip>
11
12#include "CoinTime.hpp"
13#include "CoinFileIO.hpp"
14#include "CbcModel.hpp"
15#include "CbcHeuristic.hpp"
16#include "CbcCutGenerator.hpp"
17#include "CbcBranchCut.hpp"
18#include "CglProbing.hpp"
19#include "OsiClpSolverInterface.hpp"
20#include "ClpFactorization.hpp"
21#include "OsiRowCutDebugger.hpp"
22//#############################################################################
23
24#ifdef NDEBUG
25#undef NDEBUG
26#endif
27
28//#############################################################################
29
30// Display message on stdout and stderr
31void testingMessage( const char * const msg )
32{
33    std::cout << msg;
34    //cout <<endl <<"*****************************************"
35    //     <<endl <<msg <<endl;
36}
37
38//#############################################################################
39
40static inline bool CbcTestFile(const std::string name)
41{
42    FILE *fp = fopen(name.c_str(), "r");
43    if (fp) {
44        fclose(fp);
45        return true;
46    }
47    return false;
48}
49
50//#############################################################################
51
52bool CbcTestMpsFile(std::string& fname)
53{
54    if (CbcTestFile(fname)) {
55        return true;
56    }
57    if (CbcTestFile(fname + ".mps")) {
58        fname += ".mps";
59        return true;
60    }
61    if (CbcTestFile(fname + ".MPS")) {
62        fname += ".MPS";
63        return true;
64    }
65    if (CoinFileInput::haveGzipSupport()) {
66      if (CbcTestFile(fname + ".gz")) {
67        return true;
68      }
69      if (CbcTestFile(fname + ".mps.gz")) {
70        fname += ".mps";
71        return true;
72      }
73      if (CbcTestFile(fname + ".MPS.gz")) {
74        fname += ".MPS";
75        return true;
76      }
77      if (CbcTestFile(fname + ".MPS.GZ")) {
78        fname += ".MPS";
79        return true;
80      }
81    }
82    if (CoinFileInput::haveBzip2Support()) {
83      if (CbcTestFile(fname + ".bz2")) {
84        return true;
85      }
86      if (CbcTestFile(fname + ".mps.bz2")) {
87        fname += ".mps";
88        return true;
89      }
90      if (CbcTestFile(fname + ".MPS.bz2")) {
91        fname += ".MPS";
92        return true;
93      }
94      if (CbcTestFile(fname + ".MPS.BZ2")) {
95        fname += ".MPS";
96        return true;
97      }
98    }
99    return false;
100}
101//#############################################################################
102/*
103  jjf: testSwitch -2 unitTest, -1 normal (==2)
104
105  MiplibTest might be more appropriate.
106
107  TestSwitch and stuff[6] together control how much of miplib is executed:
108    For testSwitch set to:
109      -2: solve p0033 and p0201 only (the unit test)
110      -1: solve miplib sets #0 and #1
111       0: solve nothing
112       k: execute sets j:k, where j is determined by the value of stuff[6]
113  The last parameter of PUSH_MPS specifies the test set membership.
114
115  For -miplib, -extra2 sets testSwitch, -extra3 sets stuff[6]. The command
116    cbc -extra2 -2 -miplib
117  will execute the unit test on the miplib directory.
118
119  dirMiplib should end in the directory separator character for the platform.
120
121  If you want to activate the row cut debugger for a given problem, change the
122  last parameter of the PUSH_MPS macro for the problem to true.
123
124  Returns 0 if all goes well, -1 if the Miplib directory is missing, otherwise
125     100*(number with bad objective)+(number that exceeded node limit)
126*/
127int CbcClpUnitTest (const CbcModel &saveModel, const std::string &dirMiplib,
128                    int testSwitch, const double *stuff)
129{
130  // Stop Windows popup
131  WindowsErrorPopupBlocker() ;
132  unsigned int m ;
133
134  // Do an existence check.
135  std::string test1 = dirMiplib+"p0033";
136  bool doTest = CbcTestMpsFile(test1);
137  if (!doTest) {
138    std::cout
139      << "Not doing miplib run as can't find mps files." << std::endl
140      << "Perhaps you're trying to read gzipped (.gz) files without libz?"
141      << std::endl ;
142    return (0) ;
143  }
144  int dfltPrecision = static_cast<int>(std::cout.precision()) ;
145/*
146  Set the range of problems to be tested. testSwitch = -2 is special and is
147  picked up below.
148*/
149    int loSet = 0 ;
150    int hiSet = 0 ;
151    if (testSwitch == -1) {
152      loSet = 0 ;
153      hiSet = 1 ;
154    } else if (testSwitch >= 0) {
155      loSet = static_cast<int>(stuff[6]) ;
156      hiSet = testSwitch ;
157      std::cout
158        << "Solving miplib problems in sets " << loSet
159        << ":" << hiSet << "." << std::endl ;
160    }
161/*
162  Vectors to hold test problem names and characteristics.
163*/
164    std::vector<std::string> mpsName ;
165    std::vector<int> nRows ;
166    std::vector<int> nCols ;
167    std::vector<double> objValueC ;
168    std::vector<double> objValue ;
169    std::vector<int> testSet ;
170    std::vector<int> rowCutDebugger ;
171/*
172  A macro to make the vector creation marginally readable. Parameters are
173  name, rows, columns, integer objective, continuous objective, set ID,
174  row cut debugger
175
176  To enable the row cut debugger for a given problem, change the last
177  parameter to true. Don't forget to turn it off before committing changes!
178*/
179#define PUSH_MPS(zz_mpsName_zz,\
180                 zz_nRows_zz,zz_nCols_zz,zz_objValue_zz,zz_objValueC_zz, \
181                 zz_testSet_zz, zz_rcDbg_zz) \
182  mpsName.push_back(zz_mpsName_zz) ; \
183  nRows.push_back(zz_nRows_zz) ; \
184  nCols.push_back(zz_nCols_zz) ; \
185  objValueC.push_back(zz_objValueC_zz) ; \
186  testSet.push_back(zz_testSet_zz) ; \
187  objValue.push_back(zz_objValue_zz) ; \
188  rowCutDebugger.push_back(zz_rcDbg_zz) ;
189/*
190  Push the miplib problems. Except for -2 (unitTest), push all, even if we're
191  not going to do all of them.
192*/
193  if (testSwitch == -2) {
194      PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0, false);
195      PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0, false);
196      // PUSH_MPS("flugpl", 18, 18, 1201500, 1167185.7, 0, false);
197  } else {
198/*
199  Load up the problem vector. Note that the row counts here include the
200  objective function.
201*/
202#if 1
203    PUSH_MPS("10teams", 230, 2025, 924, 917, 1, false);
204    PUSH_MPS("air03", 124, 10757, 340160, 338864.25, 0, false);
205    PUSH_MPS("air04", 823, 8904, 56137, 55535.436, 2, false);
206    PUSH_MPS("air05", 426, 7195, 26374, 25877.609, 2, false);
207    PUSH_MPS("arki001", 1048, 1388, 7580813.0459, 7579599.80787, 7, false);
208    PUSH_MPS("bell3a", 123, 133, 878430.32, 862578.64, 0, false);
209    PUSH_MPS("bell5", 91, 104, 8966406.49, 8608417.95, 1, false);
210    PUSH_MPS("blend2", 274, 353, 7.598985, 6.9156751140, 0, false);
211    PUSH_MPS("cap6000", 2176, 6000, -2451377, -2451537.325, 1, false);
212    PUSH_MPS("dano3mip", 3202, 13873, 728.1111, 576.23162474, 7, false);
213    PUSH_MPS("danoint", 664, 521, 65.67, 62.637280418, 6, false);
214    PUSH_MPS("dcmulti", 290, 548, 188182, 183975.5397, 0, false);
215    PUSH_MPS("dsbmip", 1182, 1886, -305.19817501, -305.19817501, 0, false);
216    PUSH_MPS("egout", 98, 141, 568.101, 149.589, 0, false);
217    PUSH_MPS("enigma", 21, 100, 0.0, 0.0, 0, false);
218    PUSH_MPS("fast0507", 507, 63009, 174, 172.14556668, 5, false);
219    PUSH_MPS("fiber", 363, 1298, 405935.18000, 156082.51759, 0, false);
220    PUSH_MPS("fixnet6", 478, 878, 3983, 1200.88, 1, false);
221    PUSH_MPS("flugpl", 18, 18, 1201500, 1167185.7, 0, false);
222    PUSH_MPS("gen", 780, 870, 112313, 112130.0, 0, false);
223    PUSH_MPS("gesa2", 1392, 1224, 25779856.372, 25476489.678, 1, false);
224    PUSH_MPS("gesa2_o", 1248, 1224, 25779856.372, 25476489.678, 1, false);
225    PUSH_MPS("gesa3", 1368, 1152, 27991042.648, 27833632.451, 0, false);
226    PUSH_MPS("gesa3_o", 1224, 1152, 27991042.648, 27833632.451, 0, false);
227    PUSH_MPS("gt2", 29, 188, 21166.000, 13460.233074, 0, false);
228    PUSH_MPS("harp2", 112, 2993, -73899798.00, -74353341.502, 6, false);
229    PUSH_MPS("khb05250", 101, 1350, 106940226, 95919464.0, 0, false);
230    PUSH_MPS("l152lav", 97, 1989, 4722, 4656.36, 1, false);
231    PUSH_MPS("lseu", 28, 89, 1120, 834.68, 0, false);
232    PUSH_MPS("mas74", 13, 151, 11801.18573, 10482.79528, 3, false);
233    PUSH_MPS("mas76", 12, 151, 40005.05414, 38893.9036, 2, false);
234    PUSH_MPS("misc03", 96, 160, 3360, 1910., 0, false);
235    PUSH_MPS("misc06", 820, 1808, 12850.8607, 12841.6, 0, false);
236    PUSH_MPS("misc07", 212, 260, 2810, 1415.0, 1, false);
237    PUSH_MPS("mitre", 2054, 10724, 115155, 114740.5184, 1, false);
238    PUSH_MPS("mkc", 3411, 5325, -553.75, -611.85, 7, false); // suboptimal
239    PUSH_MPS("mod008", 6, 319, 307, 290.9, 0, false);
240    PUSH_MPS("mod010", 146, 2655, 6548, 6532.08, 0, false);
241    PUSH_MPS("mod011", 4480, 10958, -54558535, -62121982.55, 2, false);
242    PUSH_MPS("modglob", 291, 422, 20740508, 20430947., 2, false);
243    PUSH_MPS("noswot", 182, 128, -43, -43.0, 6, false);
244    PUSH_MPS("nw04", 36, 87482, 16862, 16310.66667, 1, false);
245    PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0, false);
246    PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0, false);
247    PUSH_MPS("p0282", 241, 282, 258411, 176867.50, 0, false);
248    PUSH_MPS("p0548", 176, 548, 8691, 315.29, 0, false);
249    PUSH_MPS("p2756", 755, 2756, 3124, 2688.75, 0, false);
250    PUSH_MPS("pk1", 45, 86, 11.0, 0.0, 2, false);
251    PUSH_MPS("pp08a", 136, 240, 7350.0, 2748.3452381, 1, false);
252    PUSH_MPS("pp08aCUTS", 246, 240, 7350.0, 5480.6061563, 1, false);
253    PUSH_MPS("qiu", 1192, 840, -132.873137, -931.638857, 3, false);
254    PUSH_MPS("qnet1", 503, 1541, 16029.692681, 14274.102667, 0, false);
255    PUSH_MPS("qnet1_o", 456, 1541, 16029.692681, 12095.571667, 0, false);
256    PUSH_MPS("rentacar", 6803, 9557, 30356761, 28806137.644, 0, false);
257    PUSH_MPS("rgn", 24, 180, 82.1999, 48.7999, 0, false);
258    PUSH_MPS("rout", 291, 556, 1077.56, 981.86428571, 3, false);
259    PUSH_MPS("set1ch", 492, 712, 54537.75, 32007.73, 5, false);
260    PUSH_MPS("seymour", 4944, 1372, 423, 403.84647413, 7, false);
261    PUSH_MPS("seymour_1", 4944, 1372, 410.76370, 403.84647413, 5, false);
262    PUSH_MPS("stein27", 118, 27, 18, 13.0, 0, false);
263    PUSH_MPS("stein45", 331, 45, 30, 22.0, 1, false);
264    PUSH_MPS("swath", 884, 6805, 497.603, 334.4968581, 7, false);
265    PUSH_MPS("vpm1", 234, 378, 20, 15.4167, 0, false);
266    PUSH_MPS("vpm2", 234, 378, 13.75, 9.8892645972, 0, false);
267#endif
268  }
269#undef PUSH_MPS
270
271/*
272  Normally the problems are executed in order. Define RANDOM_ORDER below to
273  randomize.
274
275  #define RANDOM_ORDER
276*/
277  int which[100];
278  int nLoop = static_cast<int>(mpsName.size());
279  assert (nLoop <= 100);
280  for (int i = 0; i < nLoop; i++) which[i] = i;
281
282# ifdef RANDOM_ORDER
283  unsigned int iTime = static_cast<unsigned int>(CoinGetTimeOfDay()-1.256e9);
284  std::cout << "Time (seed) " << iTime << "." << std::endl ;
285  double sort[100];
286  CoinDrand48(true,iTime);
287  for (int i = 0; i < nLoop; i++) sort[i] = CoinDrand48();
288  CoinSort_2(sort,sort+nLoop,which);
289# endif
290
291  int problemCnt = 0;
292  for (m = 0 ; m < mpsName.size() ; m++) {
293    int setID = testSet[m];
294    if (loSet <= setID && setID <= hiSet) problemCnt++;
295  }
296
297  int numberFailures = 0;
298  int numberAttempts = 0;
299  int numProbSolved = 0;
300  double timeTaken = 0.0;
301
302//#define CLP_FACTORIZATION_INSTRUMENT
303# ifdef CLP_FACTORIZATION_INSTRUMENT
304  double timeTakenFac = 0.0;
305# endif
306
307/*
308  Open the main loop to step through the MPS problems.
309*/
310  for (unsigned int mw = 0 ; mw < mpsName.size() ; mw++) {
311    m = which[mw];
312    int setID = testSet[m];
313    // Skip if problem is not in specified problem set(s)
314    if (!(loSet <= setID && setID <= hiSet)) continue ;
315   
316    numberAttempts++;
317    std::cout << "  processing mps file: " << mpsName[m]
318              << " (" << numberAttempts << " out of "
319              << problemCnt << ")" << std::endl ;
320/*
321  Stage 1: Read the MPS and make sure the size of the constraint matrix
322           is correct.
323*/
324    CbcModel *model = new CbcModel(saveModel) ;
325
326    std::string fn = dirMiplib+mpsName[m] ;
327    if (!CbcTestMpsFile(fn)) {
328      std::cout << "ERROR: Cannot find MPS file " << fn << "." << std::endl ;
329      continue;
330    }
331    model->solver()->readMps(fn.c_str(), "") ;
332    assert(model->getNumRows() == nRows[m]) ;
333    assert(model->getNumCols() == nCols[m]) ;
334
335    // Careful! We're initialising for the benefit of other code.
336    CoinDrand48(true,1234567);
337    //printf("RAND1 %g %g\n",CoinDrand48(true,1234567),model->randomNumberGenerator()->randomDouble());
338    //printf("RAND1 %g\n",CoinDrand48(true,1234567));
339
340    // Higher limits for the serious problems.
341    int testMaximumNodes = 200000;
342    if (hiSet > 1)
343        testMaximumNodes = 20000000;
344    if (model->getMaximumNodes() > testMaximumNodes) {
345        model->setMaximumNodes(testMaximumNodes);
346    }
347/*
348  Stage 2: Call solver to solve the problem.
349*/
350
351#   ifdef CLP_FACTORIZATION_INSTRUMENT
352    extern double factorization_instrument(int type);
353    double facTime1 = factorization_instrument(0);
354    std::cout
355      << "Factorization - initial solve " << facTime1 << " seconds."
356      << std::endl ;
357      timeTakenFac += facTime1;
358#   endif
359
360    double startTime = CoinCpuTime()+CoinCpuTimeJustChildren();
361
362    // Setup specific to clp
363    OsiClpSolverInterface *siClp =
364        dynamic_cast<OsiClpSolverInterface *>(model->solver()) ;
365    ClpSimplex *modelC = NULL;
366    if (siClp) {
367      modelC = siClp->getModelPtr();
368      ClpMatrixBase * matrix = modelC->clpMatrix();
369      ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
370      if (stuff && stuff[9] && clpMatrix) {
371        // vector matrix!
372        clpMatrix->makeSpecialColumnCopy();
373      }
374
375#     ifdef JJF_ZERO
376      if (clpMatrix) {
377        int numberRows = clpMatrix->getNumRows();
378        int numberColumns = clpMatrix->getNumCols();
379        double * elements = clpMatrix->getMutableElements();
380        const int * row = clpMatrix->getIndices();
381        const CoinBigIndex * columnStart = clpMatrix->getVectorStarts();
382        const int * columnLength = clpMatrix->getVectorLengths();
383        double * smallest = new double [numberRows];
384        double * largest = new double [numberRows];
385        char * flag = new char [numberRows];
386        CoinZeroN(flag, numberRows);
387        for (int i = 0; i < numberRows; i++) {
388            smallest[i] = COIN_DBL_MAX;
389            largest[i] = 0.0;
390        }
391        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
392            bool isInteger = modelC->isInteger(iColumn);
393            CoinBigIndex j;
394            for (j = columnStart[iColumn];
395                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
396                int iRow = row[j];
397                double value = fabs(elements[j]);
398                if (!isInteger)
399                    flag[iRow] = 1;
400                smallest[iRow] = CoinMin(smallest[iRow], value);
401                largest[iRow] = CoinMax(largest[iRow], value);
402            }
403        }
404        double * rowLower = modelC->rowLower();
405        double * rowUpper = modelC->rowUpper();
406        bool changed = false;
407        for (int i = 0; i < numberRows; i++) {
408            if (flag[i] && smallest[i] > 10.0 && false) {
409                smallest[i] = 1.0 / smallest[i];
410                if (rowLower[i] > -1.0e20)
411                    rowLower[i] *= smallest[i];
412                if (rowUpper[i] < 1.0e20)
413                    rowUpper[i] *= smallest[i];
414                changed = true;
415            } else {
416                smallest[i] = 0.0;
417            }
418        }
419        if (changed) {
420            printf("SCALED\n");
421            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
422                CoinBigIndex j;
423                for (j = columnStart[iColumn];
424                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
425                    int iRow = row[j];
426                    if (smallest[iRow])
427                        elements[j] *= smallest[iRow];
428                }
429            }
430        }
431        delete [] smallest;
432        delete [] largest;
433        delete [] flag;
434      }
435#     endif    // JJF_ZERO
436
437      model->checkModel();
438      modelC->tightenPrimalBounds(0.0, 0, true);
439      model->initialSolve();
440      if (modelC->dualBound() == 1.0e10) {
441        // user did not set - so modify
442        // get largest scaled away from bound
443        ClpSimplex temp = *modelC;
444        temp.dual(0, 7);
445        double largestScaled = 1.0e-12;
446        double largest = 1.0e-12;
447        int numberRows = temp.numberRows();
448        const double * rowPrimal = temp.primalRowSolution();
449        const double * rowLower = temp.rowLower();
450        const double * rowUpper = temp.rowUpper();
451        const double * rowScale = temp.rowScale();
452        int iRow;
453        for (iRow = 0; iRow < numberRows; iRow++) {
454          double value = rowPrimal[iRow];
455          double above = value - rowLower[iRow];
456          double below = rowUpper[iRow] - value;
457          if (above < 1.0e12) {
458              largest = CoinMax(largest, above);
459          }
460          if (below < 1.0e12) {
461              largest = CoinMax(largest, below);
462          }
463          if (rowScale) {
464              double multiplier = rowScale[iRow];
465              above *= multiplier;
466              below *= multiplier;
467          }
468          if (above < 1.0e12) {
469              largestScaled = CoinMax(largestScaled, above);
470          }
471          if (below < 1.0e12) {
472              largestScaled = CoinMax(largestScaled, below);
473          }
474        }
475
476        int numberColumns = temp.numberColumns();
477        const double * columnPrimal = temp.primalColumnSolution();
478        const double * columnLower = temp.columnLower();
479        const double * columnUpper = temp.columnUpper();
480        const double * columnScale = temp.columnScale();
481        int iColumn;
482        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
483          double value = columnPrimal[iColumn];
484          double above = value - columnLower[iColumn];
485          double below = columnUpper[iColumn] - value;
486          if (above < 1.0e12) {
487              largest = CoinMax(largest, above);
488          }
489          if (below < 1.0e12) {
490              largest = CoinMax(largest, below);
491          }
492          if (columnScale) {
493              double multiplier = 1.0 / columnScale[iColumn];
494              above *= multiplier;
495              below *= multiplier;
496          }
497          if (above < 1.0e12) {
498              largestScaled = CoinMax(largestScaled, above);
499          }
500          if (below < 1.0e12) {
501              largestScaled = CoinMax(largestScaled, below);
502          }
503        }
504        std::cout << "Largest (scaled) away from bound " << largestScaled
505                  << " unscaled " << largest << std::endl;
506#       ifdef JJF_ZERO
507        modelC->setDualBound(CoinMax(1.0001e8,
508                                 CoinMin(1000.0*largestScaled,1.00001e10)));
509#       else
510        modelC->setDualBound(CoinMax(1.0001e9,
511                                 CoinMin(1000.0*largestScaled,1.0001e10)));
512#       endif
513      }
514    }    // end clp-specific setup
515/*
516  Cut passes: For small models (n < 500) always do 100 passes, if possible
517  (-100). For larger models, use minimum drop to stop (100, 20).
518*/
519    model->setMinimumDrop(CoinMin(5.0e-2,
520                      fabs(model->getMinimizationObjValue())*1.0e-3+1.0e-4));
521    if (CoinAbs(model->getMaximumCutPassesAtRoot()) <= 100) {
522        if (model->getNumCols() < 500) {
523          model->setMaximumCutPassesAtRoot(-100);
524        } else if (model->getNumCols() < 5000) {
525          model->setMaximumCutPassesAtRoot(100);
526        } else {
527          model->setMaximumCutPassesAtRoot(20);
528        }
529    }
530    // If defaults then increase trust for small models
531    if (model->numberStrong() == 5 && model->numberBeforeTrust() == 10) {
532        int numberColumns = model->getNumCols();
533        if (numberColumns <= 50) {
534            model->setNumberBeforeTrust(1000);
535        } else if (numberColumns <= 100) {
536            model->setNumberBeforeTrust(100);
537        } else if (numberColumns <= 300) {
538            model->setNumberBeforeTrust(50);
539        }
540    }
541    //if (model->getNumCols()>=500) {
542    // switch off Clp stuff
543    //model->setFastNodeDepth(-1);
544    //}
545/*
546  Activate the row cut debugger, if requested.
547*/
548    if (rowCutDebugger[m] == true) {
549      std::string probName ;
550      model->solver()->getStrParam(OsiProbName,probName) ;
551      model->solver()->activateRowCutDebugger(probName.c_str()) ;
552      if (model->solver()->getRowCutDebugger())
553        std::cout << "Row cut debugger activated for " ;
554      else
555        std::cout << "Failed to activate row cut debugger for " ;
556      std::cout << mpsName[m] << "." << std::endl ;
557    }
558    setCutAndHeuristicOptions(*model) ;
559/*
560  More clp-specific setup.
561*/
562    if (siClp) {
563#     ifdef CLP_MULTIPLE_FACTORIZATIONS
564      if (!modelC->factorization()->isDenseOrSmall()) {
565          int denseCode = stuff ? static_cast<int> (stuff[4]) : -1;
566          int smallCode = stuff ? static_cast<int> (stuff[10]) : -1;
567          if (stuff && stuff[8] >= 1) {
568              if (denseCode < 0)
569                  denseCode = 40;
570              if (smallCode < 0)
571                  smallCode = 40;
572          }
573          if (denseCode > 0)
574              modelC->factorization()->setGoDenseThreshold(denseCode);
575          if (smallCode > 0)
576              modelC->factorization()->setGoSmallThreshold(smallCode);
577          if (denseCode >= modelC->numberRows()) {
578              //printf("problem going dense\n");
579              //modelC->factorization()->goDenseOrSmall(modelC->numberRows());
580          }
581      }
582#     endif
583      if (stuff && stuff[8] >= 1) {
584        printf("Fast node size Columns %d rows %d - depth %d\n",
585               modelC->numberColumns(),modelC->numberRows(),
586               model->fastNodeDepth());
587          if (modelC->numberColumns() + modelC->numberRows() <= 10000 &&
588                  model->fastNodeDepth() == -1)
589            model->setFastNodeDepth(-10/*-9*/);
590      }
591    }
592#ifdef CONFLICT_CUTS
593    {
594      model->setCutoffAsConstraint(true); // very slow on bell5 ??
595      int moreOptions=model->moreSpecialOptions();
596      model->setMoreSpecialOptions(moreOptions|4194304);
597    }
598#endif
599/*
600  Finally, the actual call to solve the MIP with branch-and-cut.
601*/
602    model->branchAndBound();
603
604#   ifdef CLP_FACTORIZATION_INSTRUMENT
605    double facTime = factorization_instrument(0);
606    std::cout << "Factorization " << facTime << " seconds." << std::endl ,
607    timeTakenFac += facTime;
608#   endif
609
610/*
611  Stage 3: Do the statistics and check the answer.
612*/
613    double timeOfSolution = CoinCpuTime()+CoinCpuTimeJustChildren()-startTime;
614    std::cout
615      << "Cuts at root node changed objective from "
616      << model->getContinuousObjective() << " to "
617      << model->rootObjectiveAfterCuts() << std::endl ;
618    int numberGenerators = model->numberCutGenerators();
619    for (int iGenerator = 0 ; iGenerator < numberGenerators ; iGenerator++) {
620      CbcCutGenerator *generator = model->cutGenerator(iGenerator);
621#     ifdef CLIQUE_ANALYSIS
622#     ifndef CLP_INVESTIGATE
623      CglImplication *implication =
624            dynamic_cast<CglImplication*>(generator->generator());
625      if (implication) continue;
626#     endif
627#     endif
628      std::cout
629        << generator->cutGeneratorName() << " was tried "
630        << generator->numberTimesEntered() << " times and created "
631        << generator->numberCutsInTotal() << " cuts of which "
632        << generator->numberCutsActive()
633        << " were active after adding rounds of cuts";
634      if (generator->timing())
635        std::cout << " (" << generator->timeInCutGenerator() << " seconds)" ;
636      std::cout << "." << std::endl;
637    }
638    std::cout
639      << model->getNumberHeuristicSolutions()
640      << " solutions found by heuristics." << std::endl ;
641    int numberHeuristics = model->numberHeuristics();
642    for (int iHeuristic = 0 ; iHeuristic < numberHeuristics ; iHeuristic++) {
643      CbcHeuristic *heuristic = model->heuristic(iHeuristic);
644      if (heuristic->numRuns()) {
645        std::cout
646          << heuristic->heuristicName() << " was tried "
647          << heuristic->numRuns() << " times out of "
648          << heuristic->numCouldRun() << " and created "
649          << heuristic->numberSolutionsFound() << " solutions." << std::endl ;
650        }
651    }
652/*
653  Check for the correct answer.
654*/
655    if (!model->status()) {
656
657      double objActual = model->getObjValue() ;
658      double objExpect = objValue[m] ;
659      double tolerance = CoinMin(fabs(objActual),fabs(objExpect)) ;
660      tolerance = CoinMax(1.0e-5,1.0e-5*tolerance) ;
661      //CoinRelFltEq eq(1.0e-3) ;
662
663      std::cout
664        << "cbc_clp (" << mpsName[m] << ") "
665        << std::setprecision(10) << objActual ;
666      if (fabs(objActual-objExpect) < tolerance) {
667        std::cout << std::setprecision(dfltPrecision) << "; okay" ;
668          numProbSolved++;
669      } else  {
670        std::cout
671          << " != " << objExpect << std::setprecision(dfltPrecision)
672          << "; error = " << fabs(objExpect-objActual) ;
673          numberFailures++;
674        //#ifdef COIN_DEVELOP
675        //abort();
676        //#endif
677      }
678    } else {
679      std::cout
680        << "cbc_clp (" << mpsName[m] << ") status not optimal; "
681        << "assuming too many nodes" ;
682    }
683    timeTaken += timeOfSolution;
684    std::cout
685      << " -- (" << model->getNodeCount() << " n / "
686      << model->getIterationCount() << " i / "
687      << timeOfSolution << " s) (subtotal " << timeTaken << " seconds)"
688      << std::endl;
689    delete model;
690  }
691/*
692  End main loop on MPS problems. Print a summary and calculate the return
693  value.
694*/
695  int returnCode = 0;
696  std::cout
697    << "cbc_clp solved " << numProbSolved << " out of " << numberAttempts;
698  int numberOnNodes = numberAttempts-numProbSolved-numberFailures;
699  if (numberFailures || numberOnNodes) {
700    if (numberOnNodes) {
701      std::cout << " (" << numberOnNodes << " stopped on nodes)";
702      returnCode = numberOnNodes;
703    }
704    if (numberFailures) {
705      std::cout << " (" << numberFailures << " gave bad answer!)";
706      returnCode += 100*numberFailures;
707    }
708  }
709  std::cout
710    << " and took " << timeTaken << " seconds." << std::endl;
711
712  if (testSwitch == -2) {
713    if (numberFailures || numberOnNodes) {
714      std::cout << "****** Unit Test failed." << std::endl ;
715      std::cerr << "****** Unit Test failed." << std::endl ;
716    } else {
717      std::cerr << "****** Unit Test succeeded." << std::endl ;
718    }
719  }
720# ifdef CLP_FACTORIZATION_INSTRUMENT
721  std::cout
722    << "Total factorization time " << timeTakenFac << "seconds." << std::endl ;
723# endif
724  return (returnCode) ;
725}
726
Note: See TracBrowser for help on using the repository browser.