Changeset 1622 for trunk/Cbc


Ignore:
Timestamp:
Mar 25, 2011 1:20:10 PM (9 years ago)
Author:
lou
Message:

Clean up CbcClpUnitTest?. Easy activation for row cut debugger. Documentation.

Location:
trunk/Cbc/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcSolver.cpp

    r1621 r1622  
    10991099
    11001100int CbcClpUnitTest (const CbcModel & saveModel,
    1101                     std::string& dirMiplib, int testSwitch,
    1102                     double * stuff);
     1101                    const std::string& dirMiplib, int testSwitch,
     1102                    const double * stuff);
    11031103
    11041104int CbcMain1 (int argc, const char *argv[],
  • trunk/Cbc/src/unitTestClp.cpp

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