Changeset 764 for branches


Ignore:
Timestamp:
Aug 19, 2007 5:23:49 PM (12 years ago)
Author:
ladanyi
Message:

committing merging w/ trunk of Clp, Osi, Cbc

Location:
branches/BSP/trunk/Cbc
Files:
12 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/BSP/trunk/Cbc/src/CbcEventHandler.hpp

    r706 r764  
    8888    treeStatus,
    8989    /*! A solution has been found. */
    90     solution
     90    solution,
     91    /*! A heuristic solution has been found. */
     92    heuristicSolution
    9193  } ;
    9294
  • branches/BSP/trunk/Cbc/src/CbcHeuristic.cpp

    r732 r764  
    141141    CglPreProcess process;
    142142    /* Do not try and produce equality cliques and
    143        do up to 5 passes */
     143       do up to 2 passes */
    144144    if (logLevel<=1)
    145145      process.messageHandler()->setLogLevel(0);
    146     OsiSolverInterface * solver2= process.preProcess(*solver);
     146    OsiSolverInterface * solver2= process.preProcessNonDefault(*solver,false,2);
    147147    if (!solver2) {
    148148      if (logLevel>1)
     
    153153      double before = solver->getNumRows()+solver->getNumCols();
    154154      double after = solver2->getNumRows()+solver2->getNumCols();
    155       if (logLevel>1)
    156         printf("before %d rows %d columns, after %d rows %d columns\n",
    157                solver->getNumRows(),solver->getNumCols(),
    158                solver2->getNumRows(),solver2->getNumCols());
    159       if (after>fractionSmall_*before)
    160         return 0;
     155      char generalPrint[100];
     156      if (after>fractionSmall_*before&&after>300) {
     157        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
     158                solver->getNumRows(),solver->getNumCols(),
     159                solver2->getNumRows(),solver2->getNumCols());
     160        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     161          << generalPrint
     162          <<CoinMessageEol;
     163        return -1;
     164      } else {
     165        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns",
     166                solver->getNumRows(),solver->getNumCols(),
     167                solver2->getNumRows(),solver2->getNumCols());
     168        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     169          << generalPrint
     170          <<CoinMessageEol;
     171      }
    161172      solver2->resolve();
    162173      CbcModel model(*solver2);
     
    191202      model.setMaximumCutPassesAtRoot(CoinMin(20,model_->getMaximumCutPassesAtRoot()));
    192203      model.setParentModel(*model_);
     204      model.setOriginalColumns(process.originalColumns());
    193205      model.branchAndBound();
    194206      if (logLevel>1)
  • branches/BSP/trunk/Cbc/src/CbcHeuristicFPump.cpp

    r723 r764  
    1515#include "CbcBranchActual.hpp"
    1616#include "CoinHelperFunctions.hpp"
     17#include "CoinWarmStartBasis.hpp"
    1718#include "CoinTime.hpp"
     19#include "CbcEventHandler.hpp"
    1820
    1921
     
    3234   maximumRetries_(1),
    3335   accumulate_(0),
     36   fixOnReducedCosts_(1),
    3437   roundExpensive_(false)
    3538{
     
    5255   maximumRetries_(1),
    5356   accumulate_(0),
     57   fixOnReducedCosts_(1),
    5458   roundExpensive_(roundExpensive)
    5559{
     
    8892  else
    8993    fprintf(fp,"4  heuristicFPump.setAccumulate(%d);\n",accumulate_);
     94  if (fixOnReducedCosts_!=other.fixOnReducedCosts_)
     95    fprintf(fp,"3  heuristicFPump.setFixOnReducedCosts(%d);\n",fixOnReducedCosts_);
     96  else
     97    fprintf(fp,"4  heuristicFPump.setFixOnReducedCosts(%d);\n",fixOnReducedCosts_);
    9098  if (maximumTime_!=other.maximumTime_)
    9199    fprintf(fp,"3  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
     
    134142  maximumRetries_(rhs.maximumRetries_),
    135143  accumulate_(rhs.accumulate_),
     144  fixOnReducedCosts_(rhs.fixOnReducedCosts_),
    136145  roundExpensive_(rhs.roundExpensive_)
    137146{
     
    155164    maximumRetries_ = rhs.maximumRetries_;
    156165    accumulate_ = rhs.accumulate_;
     166    fixOnReducedCosts_ = rhs.fixOnReducedCosts_;
    157167    roundExpensive_ = rhs.roundExpensive_;
    158168  }
     
    187197    return 0;
    188198  // probably a good idea
    189   if (model_->getSolutionCount()) return 0;
     199  //if (model_->getSolutionCount()) return 0;
    190200  // loop round doing repeated pumps
    191201  double cutoff;
     
    253263  }
    254264  double time1 = CoinCpuTime();
     265  model_->solver()->resolve();
     266  if (cutoff<1.0e50&&false) {
     267    // Fix on djs
     268    double direction = model_->solver()->getObjSense() ;
     269    double gap = cutoff - model_->solver()->getObjValue()*direction ;
     270    double tolerance;
     271    model_->solver()->getDblParam(OsiDualTolerance,tolerance) ;
     272    if (gap>0.0) {
     273      gap += 100.0*tolerance;
     274      int nFix=model_->solver()->reducedCostFix(gap);
     275      printf("dj fixing fixed %d variables\n",nFix);
     276    }
     277  }
     278  CoinWarmStartBasis saveBasis;
     279  CoinWarmStartBasis * basis =
     280    dynamic_cast<CoinWarmStartBasis *>(model_->solver()->getWarmStart()) ;
     281  if (basis) {
     282    saveBasis = * basis;
     283    delete basis;
     284  }
     285  double continuousObjectiveValue = model_->solver()->getObjValue()*model_->solver()->getObjSense();
     286  double * firstPerturbedObjective = NULL;
     287  double * firstPerturbedSolution = NULL;
     288  assert (model_->solver()->isProvenOptimal());
    255289  if (when_>=11&&when_<=15) {
    256290    fixInternal = when_ >11&&when_<15;
     
    264298    usedColumn = new char [numberColumns];
    265299    memset(usedColumn,0,numberColumns);
    266     model_->solver()->resolve();
    267     assert (model_->solver()->isProvenOptimal());
    268300    lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(),numberColumns);
    269301  }
     
    271303  int totalNumberPasses=0;
    272304  int numberTries=0;
    273   while (true) {
     305  CoinWarmStartBasis bestBasis;
     306  bool exitAll=false;
     307  double saveBestObjective = model_->getMinimizationObjValue();
     308  while (!exitAll) {
    274309    int numberPasses=0;
    275310    numberTries++;
    276311    // Clone solver - otherwise annoys root node computations
    277312    OsiSolverInterface * solver = model_->solver()->clone();
     313    if (CoinMin(fakeCutoff_,cutoff)<1.0e50) {
     314      // Fix on djs
     315      double direction = solver->getObjSense() ;
     316      double gap = CoinMin(fakeCutoff_,cutoff) - solver->getObjValue()*direction ;
     317      double tolerance;
     318      solver->getDblParam(OsiDualTolerance,tolerance) ;
     319      if (gap>0.0&&(fixOnReducedCosts_==1||(numberTries==1&&fixOnReducedCosts_==2))) {
     320        gap += 100.0*tolerance;
     321        int nFix=solver->reducedCostFix(gap);
     322        if (nFix) {
     323          sprintf(pumpPrint,"Reduced cost fixing fixed %d variables on pass %d",nFix,numberTries);
     324          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     325            << pumpPrint
     326            <<CoinMessageEol;
     327          pumpPrint[0]='\0';
     328        }
     329      }
     330    }
    278331    // if cutoff exists then add constraint
    279332    bool useCutoff = (fabs(cutoff)<1.0e20&&(fakeCutoff_!=COIN_DBL_MAX||numberTries>1));
     
    295348        }
    296349      }
    297       solver->addRow(nel,which,els,-COIN_DBL_MAX,cutoff);
     350      double offset;
     351      solver->getDblParam(OsiObjOffset,offset);
     352#ifdef COIN_DEVELOP
     353      if (offset)
     354        printf("CbcHeuristicFPump obj offset %g\n",offset);
     355#endif
     356      solver->addRow(nel,which,els,-COIN_DBL_MAX,cutoff+offset);
    298357      delete [] which;
    299358      delete [] els;
     
    356415    while (!finished) {
    357416      returnCode=0;
     417      if (model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
     418        exitAll=true;
     419        break;
     420      }
    358421      // see what changed
    359422      if (usedColumn) {
     
    370433      memcpy(newSolution,solution,numberColumns*sizeof(double));
    371434      int flip;
     435      if (numberPasses==0&&false) {
     436        // always use same seed
     437        CoinSeedRandom(987654321);
     438      }
    372439      returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
    373440                          pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
     441      if (numberPasses==0&&false) {
     442        // Make sure random will be different
     443        for (i=1;i<numberTries;i++)
     444          CoinDrand48();
     445      }
    374446      numberPasses++;
    375447      if (returnCode) {
     
    389461        newLineNeeded=false;
    390462        if (newSolutionValue<solutionValue) {
    391           double saveValue = newSolutionValue;
     463          double saveValue = solutionValue;
    392464          if (!doGeneral) {
    393465            int numberLeft=0;
     
    406478              returnCode = smallBranchAndBound(solver,numberNodes_,newSolution,newSolutionValue,
    407479                                               solutionValue,"CbcHeuristicFpump");
     480              if (returnCode<0)
     481                returnCode=0; // returned on size
    408482              if ((returnCode&2)!=0) {
    409483                // could add cut
    410484                returnCode &= ~2;
    411485              }
    412             }
    413           }
    414           if (returnCode) {
     486              if (returnCode!=1)
     487                newSolutionValue=saveValue;
     488            }
     489          }
     490          if (returnCode&&newSolutionValue<saveValue) {
    415491            memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     492            solutionFound=true;
     493            CoinWarmStartBasis * basis =
     494              dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     495            if (basis) {
     496              bestBasis = * basis;
     497              delete basis;
     498              CbcEventHandler * handler = model_->getEventHandler();
     499              if (handler) {
     500                double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
     501                double saveObjectiveValue = model_->getMinimizationObjValue();
     502                model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
     503                int action = handler->event(CbcEventHandler::heuristicSolution);
     504                if (saveOldSolution) {
     505                  model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
     506                  delete [] saveOldSolution;
     507                }
     508                if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
     509                  exitAll=true; // exit
     510                  break;
     511                }
     512              }
     513            }
    416514            if ((accumulate_&1)!=0)
    417515              model_->incrementUsed(betterSolution); // for local search
     
    426524            pumpPrint[0]='\0';
    427525          } else {
    428             sprintf(pumpPrint+strlen(pumpPrint)," - not good enough after mini branch and bound");
     526            sprintf(pumpPrint+strlen(pumpPrint)," - mini branch and bound could not fix general integers");
    429527            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
    430528              << pumpPrint
     
    458556          if (matched) break;
    459557        }
     558        int numberPerturbed=0;
    460559        if (matched || numberPasses%100 == 0) {
    461560          // perturbation
    462561          sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
    463562          newLineNeeded=true;
     563          double factorX[10]={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
     564          double factor=1.0;
     565          double target=-1.0;
     566          double * randomX = new double [numberIntegers];
     567          for (i=0;i<numberIntegers;i++)
     568            randomX[i] = max(0.0,CoinDrand48()-0.3);
     569          for (int k=0;k<10;k++) {
     570#ifdef COIN_DEVELOP_x
     571            printf("kpass %d\n",k);
     572#endif
     573            int numberX[10]={0,0,0,0,0,0,0,0,0,0};
     574            for (i=0;i<numberIntegers;i++) {
     575              int iColumn = integerVariable[i];
     576              double value = randomX[i];
     577              double difference = fabs(solution[iColumn]-newSolution[iColumn]);
     578              for (int j=0;j<10;j++) {
     579                if (difference+value*factorX[j]>0.5)
     580                  numberX[j]++;
     581              }
     582            }
     583            if (target<0.0) {
     584              if (numberX[9]<=200)
     585                break; // not very many changes
     586              target=CoinMax(200.0,CoinMin(0.05*numberX[9],1000.0));
     587            }
     588            int iX=-1;
     589            int iBand=-1;
     590            for (i=0;i<10;i++) {
     591#ifdef COIN_DEVELOP_x
     592              printf("** %d changed at %g\n",numberX[i],factorX[i]);
     593#endif
     594              if (numberX[i]>=target&&numberX[i]<2.0*target&&iX<0)
     595                iX=i;
     596              if (iBand<0&&numberX[i]>target) {
     597                iBand=i;
     598                factor=factorX[i];
     599              }
     600            }
     601            if (iX>=0) {
     602              factor=factorX[iX];
     603              break;
     604            } else {
     605              assert (iBand>=0);
     606              double hi = factor;
     607              double lo = (iBand>0) ? factorX[iBand-1] : 0.0;
     608              double diff = (hi-lo)/9.0;
     609              for (i=0;i<10;i++) {
     610                factorX[i]=lo;
     611                lo += diff;
     612              }
     613            }
     614          }
    464615          for (i=0;i<numberIntegers;i++) {
    465616            int iColumn = integerVariable[i];
    466             double value = max(0.0,CoinDrand48()-0.3);
     617            double value = randomX[i];
    467618            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
    468             if (difference+value>0.5) {
     619            if (difference+value*factor>0.5) {
     620              numberPerturbed++;
    469621              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
    470622                newSolution[iColumn] += 1.0;
     
    480632            }
    481633          }
     634          delete [] randomX;
    482635        } else {
    483636          for (j=NUMBER_OLD-1;j>0;j--) {
     
    490643        double offset=0.0;
    491644        double costValue = (1.0-scaleFactor)*solver->getObjSense();
    492        
     645        int numberChanged=0;
     646        const double * oldObjective = solver->getObjCoefficients();
    493647        for (i=0;i<numberColumns;i++) {
    494648          // below so we can keep original code and allow for objective
     
    502656            continue;
    503657          }
     658          double newValue=0.0;
    504659          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
    505             solver->setObjCoeff(iColumn,costValue+scaleFactor*saveObjective[iColumn]);
     660            newValue = costValue+scaleFactor*saveObjective[iColumn];
    506661          } else {
    507662            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
    508               solver->setObjCoeff(iColumn,-costValue+scaleFactor*saveObjective[iColumn]);
    509             } else {
    510               solver->setObjCoeff(iColumn,0.0);
    511             }
    512           }
     663              newValue = -costValue+scaleFactor*saveObjective[iColumn];
     664            }
     665          }
     666          if (newValue!=oldObjective[iColumn])
     667            numberChanged++;
     668          solver->setObjCoeff(iColumn,newValue);
    513669          offset += costValue*newSolution[iColumn];
    514670        }
     
    555711            if (newSolutionValue<solutionValue) {
    556712              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     713              CoinWarmStartBasis * basis =
     714                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     715              solutionFound=true;
     716              if (basis) {
     717                bestBasis = * basis;
     718                delete basis;
     719                CbcEventHandler * handler = model_->getEventHandler();
     720                if (handler) {
     721                  double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
     722                  double saveObjectiveValue = model_->getMinimizationObjValue();
     723                  model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
     724                  int action = handler->event(CbcEventHandler::heuristicSolution);
     725                  if (saveOldSolution) {
     726                    model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
     727                    delete [] saveOldSolution;
     728                  }
     729                  if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
     730                    exitAll=true; // exit
     731                    break;
     732                  }
     733                }
     734              }
    557735              if ((accumulate_&1)!=0)
    558736                model_->incrementUsed(betterSolution); // for local search
     
    570748            solver->setWarmStart(&dummy);
    571749          }
     750#ifdef COIN_DEVELOP
     751          printf("%d perturbed out of %d columns (%d changed)\n",numberPerturbed,numberColumns,numberChanged);
     752#endif
     753          bool takeHint;
     754          OsiHintStrength strength;
     755          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
     756          if (numberPerturbed>numberColumns)
     757            solver->setHintParam(OsiDoDualInResolve,true); // dual may be better
     758          if (numberTries>1&&numberPasses==1&&false) {
     759            // use basis from first time
     760            solver->setWarmStart(&saveBasis);
     761            // and objective
     762            solver->setObjective(firstPerturbedObjective);
     763            // and solution
     764            solver->setColSolution(firstPerturbedSolution);
     765          }
    572766          solver->resolve();
     767          if (numberTries==1&&numberPasses==1&&false) {
     768            // save basis
     769            CoinWarmStartBasis * basis =
     770              dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     771            if (basis) {
     772              saveBasis = * basis;
     773              delete basis;
     774            }
     775            firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(),numberColumns);
     776            firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(),numberColumns);
     777          }
     778          solver->setHintParam(OsiDoDualInResolve,takeHint);
     779#ifdef COIN_DEVELOP
     780          {
     781            double newSolutionValue = -saveOffset;
     782            const double * newSolution = solver->getColSolution();
     783            for (  i=0 ; i<numberColumns ; i++ )
     784              newSolutionValue += saveObjective[i]*newSolution[i];
     785            newSolutionValue *= direction;
     786            printf("took %d iterations - true obj %g\n",solver->getIterationCount(),newSolutionValue);
     787          }
     788#endif
    573789          assert (solver->isProvenOptimal());
    574790          // in case very dubious solver
     
    675891      scaleFactor *= weightFactor_;
    676892    } // END WHILE
    677     if (!solutionFound) {
     893    if (!solutionFound)
    678894      sprintf(pumpPrint+strlen(pumpPrint),"No solution found this major pass");
     895    if (strlen(pumpPrint)) {
    679896      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
    680897        << pumpPrint
     
    687904    delete [] oldSolution;
    688905    delete [] saveObjective;
    689     if (usedColumn) {
     906    if (usedColumn&&!exitAll) {
    690907      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
    691908      const double * colLower = newSolver->getColLower();
     
    764981      returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
    765982                                       cutoff,"CbcHeuristicLocalAfterFPump");
     983      if (returnCode<0)
     984        returnCode=0; // returned on size - could try changing
    766985      if ((returnCode&2)!=0) {
    767986        // could add cut
    768987        returnCode &= ~2;
    769988      }
    770       if (returnCode) {
     989      if (returnCode&&newSolutionValue<saveValue) {
    771990        sprintf(pumpPrint+strlen(pumpPrint),"Mini branch and bound improved solution from %g to %g (%.2f seconds)",
    772991                saveValue,newSolutionValue,model_->getCurrentSeconds());
     
    8101029        solutionValue=newSolutionValue;
    8111030        solutionFound=true;
     1031        CoinWarmStartBasis * basis =
     1032          dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
     1033        if (basis) {
     1034          bestBasis = * basis;
     1035          delete basis;
     1036          CbcEventHandler * handler = model_->getEventHandler();
     1037          if (handler) {
     1038            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
     1039            double saveObjectiveValue = model_->getMinimizationObjValue();
     1040            model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
     1041            int action = handler->event(CbcEventHandler::heuristicSolution);
     1042            if (saveOldSolution) {
     1043              model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
     1044              delete [] saveOldSolution;
     1045            }
     1046            if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
     1047              exitAll=true; // exit
     1048              break;
     1049            }
     1050          }
     1051        }
    8121052      } else {
    8131053        sprintf(pumpPrint+strlen(pumpPrint),"Mini branch and bound did not improve solution (%.2f seconds)",
     
    8221062    if (solutionFound) finalReturnCode=1;
    8231063    cutoff = CoinMin(cutoff,solutionValue);
    824     if (numberTries>=maximumRetries_||!solutionFound) {
     1064    if (numberTries>=maximumRetries_||!solutionFound||exitAll||cutoff<continuousObjectiveValue+1.0e-7) {
    8251065      break;
    826     } else if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
     1066    } else {
    8271067      solutionFound=false;
    828       double gap = relativeIncrement_*fabs(solutionValue);
    829       cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
     1068      if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
     1069        double gap = relativeIncrement_*fabs(solutionValue);
     1070        cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
     1071      } else {
     1072        cutoff -= 0.1*(cutoff-continuousObjectiveValue);
     1073      }
     1074      if (cutoff<continuousObjectiveValue)
     1075        break;
    8301076      sprintf(pumpPrint+strlen(pumpPrint),"Round again with cutoff of %g",cutoff);
    8311077      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     
    8361082        memset(usedColumn,0,numberColumns);
    8371083      totalNumberPasses += numberPasses-1;
    838     } else {
    839       break;
    8401084    }
    8411085  }
     
    8671111    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
    8681112                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
     1113    if (returnCode<0)
     1114      returnCode=0; // returned on size
    8691115    if ((returnCode&2)!=0) {
    8701116      // could add cut
     
    8851131  delete [] closestSolution;
    8861132  delete [] integerVariable;
     1133  delete [] firstPerturbedObjective;
     1134  delete [] firstPerturbedSolution;
    8871135  sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
    8881136          model_->getCurrentSeconds(),CoinCpuTime()-time1);
     
    8901138    << pumpPrint
    8911139    <<CoinMessageEol;
     1140  if (bestBasis.getNumStructural())
     1141    model_->setBestSolutionBasis(bestBasis);
     1142  model_->setMinimizationObjValue(saveBestObjective);
    8921143  return finalReturnCode;
    8931144}
  • branches/BSP/trunk/Cbc/src/CbcHeuristicFPump.hpp

    r746 r764  
    121121  inline int accumulate() const
    122122  { return accumulate_;}
     123  /**  Set whether to fix variables on known solution
     124       0 - do not fix
     125       1 - fix integers on reduced costs
     126       2 - fix integers on reduced costs but only on entry
     127  */
     128  inline void setFixOnReducedCosts(int value)
     129  { fixOnReducedCosts_=value;}
     130  /// Get reduced cost option
     131  inline int fixOnReducedCosts() const
     132  { return fixOnReducedCosts_;}
    123133
    124134protected:
     
    156166  */
    157167  int accumulate_;
     168  /**  Set whether to fix variables on known solution
     169       0 - do not fix
     170       1 - fix integers on reduced costs
     171       2 - fix integers on reduced costs but only on entry
     172  */
     173  int fixOnReducedCosts_;
    158174  /// If true round to expensive
    159175  bool roundExpensive_;
  • branches/BSP/trunk/Cbc/src/CbcHeuristicGreedy.cpp

    r687 r764  
    731731    int returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
    732732                                         solutionValue,"CbcHeuristicGreedy");
     733    if (returnCode<0)
     734      returnCode=0; // returned on size
    733735    if ((returnCode&2)!=0) {
    734736      // could add cut
  • branches/BSP/trunk/Cbc/src/CbcHeuristicLocal.cpp

    r708 r764  
    156156  int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,objectiveValue,
    157157                                         objectiveValue,"CbcHeuristicLocal");
     158  if (returnCode<0)
     159    returnCode=0; // returned on size
    158160  if ((returnCode&2)!=0) {
    159161    // could add cut
     
    547549        // new solution
    548550        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     551        CoinWarmStartBasis * basis =
     552          dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     553        if (basis) {
     554          model_->setBestSolutionBasis(* basis);
     555          delete basis;
     556        }
    549557        returnCode=1;
    550558        solutionValue = newSolutionValue + bestChange;
  • branches/BSP/trunk/Cbc/src/CbcHeuristicRINS.cpp

    r640 r764  
    203203      returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
    204204                                         model_->getCutoff(),"CbcHeuristicRINS");
     205      if (returnCode<0)
     206        returnCode=0; // returned on size
    205207      if ((returnCode&1)!=0)
    206208        numberSuccesses_++;
  • branches/BSP/trunk/Cbc/src/CbcMain.cpp

    r742 r764  
    859859            }
    860860          } else if (valid==1) {
    861             abort();
     861            std::cout<<" is illegal for double parameter "<<parameters[iParam].name()<<" value remains "<<
     862              parameters[iParam].doubleValue()<<std::endl;
    862863          } else {
    863864            switch(type) {
     
    883884            parameters[iParam].setIntParameter(*model,value);
    884885          } else if (valid==1) {
    885             abort();
     886            std::cout<<" is illegal for integer parameter "<<parameters[iParam].name()<<" value remains "<<
     887              parameters[iParam].intValue()<<std::endl;
    886888          } else {
    887889            std::cout<<parameters[iParam].name()<<" has value "<<
  • branches/BSP/trunk/Cbc/src/CbcModel.cpp

    r736 r764  
    513513  strongInfo_[2]=0;
    514514  numberStrongIterations_ = 0;
     515  // Initialize random seed
     516  CoinSeedRandom(1234567);
    515517#ifndef NDEBUG
    516518  {
     519#ifdef COIN_DEVELOP
     520    double big = 1.0e10;
     521#else
     522    double big = 1.0e20;
     523#endif
    517524    int i;
    518525    int n = solver_->getNumCols();
     
    520527    const double *upper = solver_->getColUpper() ;
    521528    for (i=0;i<n;i++) {
    522       assert (lower[i]<1.0e10);
    523       assert (upper[i]>-1.0e10);
     529      assert (lower[i]<big);
     530      assert (upper[i]>-big);
    524531    }
    525532    n = solver_->getNumRows();
     
    527534    upper = solver_->getRowUpper() ;
    528535    for (i=0;i<n;i++) {
    529       assert (lower[i]<1.0e10);
    530       assert (upper[i]>-1.0e10);
     536      assert (lower[i]<big);
     537      assert (upper[i]>-big);
    531538    }
    532539  }
     
    691698  eventHappened_=false;
    692699  CbcEventHandler *eventHandler = getEventHandler() ;
     700  if (eventHandler)
     701    eventHandler->setModel(this);
    693702  // set up for probing
    694703  probingInfo_ = new CglTreeProbingInfo(solver_);
     
    936945*/
    937946  bestObjective_ = CoinMin(bestObjective_,1.0e50) ;
    938   numberSolutions_ = 0 ;
    939   stateOfSearch_=0;
    940   numberHeuristicSolutions_ = 0 ;
     947  if (!bestSolution_) {
     948    numberSolutions_ = 0 ;
     949    numberHeuristicSolutions_ = 0 ;
     950  }
     951  stateOfSearch_ = 0;
    941952  // Everything is minimization
    942953  {
     
    10821093  numberDJFixed_=0.0;
    10831094  // Do heuristics
    1084   {
    1085     double * newSolution = new double [numberColumns] ;
    1086     double heuristicValue = getCutoff() ;
    1087     int found = -1; // no solution found
    1088 
    1089     currentPassNumber_ = 1; // so root heuristics will run
    1090     int i;
    1091     for (i = 0;i<numberHeuristics_;i++) {
    1092       // see if heuristic will do anything
    1093       double saveValue = heuristicValue ;
    1094       int ifSol = heuristic_[i]->solution(heuristicValue,
    1095                                           newSolution);
    1096       if (ifSol>0) {
    1097         // better solution found
    1098         found = i ;
    1099         incrementUsed(newSolution);
    1100         // increment number of solutions so other heuristics can test
    1101         numberSolutions_++;
    1102         numberHeuristicSolutions_++;
    1103       } else {
    1104         heuristicValue = saveValue ;
    1105       }
    1106     }
    1107     currentPassNumber_ = 0;
    1108     /*
    1109       Did any of the heuristics turn up a new solution? Record it before we free
    1110       the vector.
    1111     */
    1112     if (found >= 0) {
    1113       // For compiler error on terra cluster!
    1114       if (found<numberHeuristics_)
    1115         lastHeuristic_ = heuristic_[found];
    1116       else
    1117         lastHeuristic_ = heuristic_[0];
    1118       setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    1119       CbcTreeLocal * tree
    1120           = dynamic_cast<CbcTreeLocal *> (tree_);
    1121       if (tree)
    1122         tree->passInSolution(bestSolution_,heuristicValue);
    1123       for (i = 0;i<numberHeuristics_;i++) {
    1124         // delete FPump
    1125         CbcHeuristicFPump * pump
    1126           = dynamic_cast<CbcHeuristicFPump *> (heuristic_[i]);
    1127         if (pump) {
    1128           delete pump;
    1129           numberHeuristics_ --;
    1130           for (int j=i;j<numberHeuristics_;j++)
    1131             heuristic_[j] = heuristic_[j+1];
    1132         }
    1133       }
    1134     }
    1135     delete [] newSolution ;
    1136   }
     1095  doHeuristicsAtRoot();
    11371096  statistics_ = NULL;
    11381097  // Do on switch
     
    11591118
    11601119      if (numberUnsatisfied)   {
    1161         feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
    1162                                  NULL);
     1120        // User event
     1121        if (!eventHappened_)
     1122          feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
     1123                                   NULL);
     1124        else
     1125          feasible=false;
    11631126      } else if (solverCharacteristics_->solutionAddsCuts()||
    11641127                 solverCharacteristics_->alwaysTryCutsAtRootNode()) {
     
    11901153    feasible = false;
    11911154  }
     1155  // User event
     1156  if (eventHappened_)
     1157    feasible=false;
    11921158/*
    11931159  We've taken the continuous relaxation as far as we can. Time to branch.
     
    22612227      pthread_cond_signal(threadInfo[i].condition2); // unlock
    22622228      //if (!stopped)
    2263         pthread_join(threadId[i],NULL);
     2229      //pthread_join(threadId[i],NULL);
     2230      int returnCode;
     2231      returnCode=pthread_join(threadId[i],NULL);
     2232      assert (!returnCode);
    22642233        //else
    22652234        //pthread_kill(threadId[i]); // kill rather than try and synchronize
     
    30883057{
    30893058  // resize best solution if exists
    3090   if (bestSolution_) {
     3059  if (bestSolution_&&solver&&solver_) {
    30913060    int nOld = solver_->getNumCols();
    30923061    int nNew = solver->getNumCols();
     
    31133082  { delete emptyWarmStart_  ;
    31143083    emptyWarmStart_ = 0 ; }
     3084  bestSolutionBasis_ = CoinWarmStartBasis();
    31153085/*
    31163086  Initialize integer variable vector.
     
    32493219  lastHeuristic_ = NULL;
    32503220  if (rhs.eventHandler_)
    3251   { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; }
     3221    { eventHandler_ = rhs.eventHandler_->clone() ; }
    32523222  else
    32533223  { eventHandler_ = NULL ; }
     
    33893359    addedCuts_ = NULL;
    33903360  }
     3361  bestSolutionBasis_ = rhs.bestSolutionBasis_;
    33913362  nextRowCut_ = NULL;
    33923363  currentNode_ = NULL;
     
    35693540      delete eventHandler_ ;
    35703541    if (rhs.eventHandler_)
    3571     { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; }
     3542      { eventHandler_ = rhs.eventHandler_->clone() ; }
    35723543    else
    35733544    { eventHandler_ = NULL ; }
     
    36503621      addedCuts_ = NULL;
    36513622    }
     3623    bestSolutionBasis_ = rhs.bestSolutionBasis_;
    36523624    nextRowCut_ = NULL;
    36533625    currentNode_ = NULL;
     
    49074879            // add to global list
    49084880            OsiRowCut newCut(*thisCut);
    4909             newCut.setGloballyValid(2);
     4881            newCut.setGloballyValid(true);
    49104882            newCut.mutableRow().setTestForDuplicateIndex(false);
    49114883            globalCuts_.insert(newCut) ;
     
    49274899              // add to global list
    49284900              OsiRowCut newCut(*thisCut);
    4929               newCut.setGloballyValid(2);
     4901              newCut.setGloballyValid(true);
    49304902              newCut.mutableRow().setTestForDuplicateIndex(false);
    49314903              globalCuts_.insert(newCut) ;
     
    49394911            // add to global list
    49404912            OsiColCut newCut(*thisCut);
    4941             newCut.setGloballyValid(2);
     4913            newCut.setGloballyValid(true);
    49424914            globalCuts_.insert(newCut) ;
    49434915          }
     
    51625134            // add to global list
    51635135            OsiRowCut newCut(*thisCut);
    5164             newCut.setGloballyValid(2);
     5136            newCut.setGloballyValid(true);
    51655137            newCut.mutableRow().setTestForDuplicateIndex(false);
    51665138            globalCuts_.insert(newCut) ;
     
    51735145            // add to global list
    51745146            OsiColCut newCut(*thisCut);
    5175             newCut.setGloballyValid(2);
     5147            newCut.setGloballyValid(true);
    51765148            globalCuts_.insert(newCut) ;
    51775149          }
     
    53475319*/
    53485320    int numberRowsNow = solver_->getNumRows() ;
     5321#ifndef NDEBUG
    53495322    assert(numberRowsNow == numberRowsAtStart+lastNumberCuts) ;
     5323#else
     5324    // ? maybe clue to threaded problems
     5325    if(numberRowsNow != numberRowsAtStart+lastNumberCuts) {
     5326      fprintf(stderr,"*** threaded error - numberRowsNow(%d) != numberRowsAtStart(%d)+lastNumberCuts(%d)\n",
     5327              numberRowsNow,numberRowsAtStart,lastNumberCuts);
     5328      fprintf(stdout,"*** threaded error - numberRowsNow(%d) != numberRowsAtStart(%d)+lastNumberCuts(%d)\n",
     5329              numberRowsNow,numberRowsAtStart,lastNumberCuts);
     5330      abort();
     5331    }
     5332#endif
    53505333    int numberToAdd = theseCuts.sizeRowCuts() ;
    53515334    numberNewCuts_ = lastNumberCuts+numberToAdd ;
     
    73277310        solver_->setWarmStart(slack);
    73287311        delete slack ;
     7312      } else {
     7313        if (bestSolutionBasis_.getNumStructural()==solver_->getNumCols()&&
     7314            bestSolutionBasis_.getNumArtificial()==solver_->getNumRows())
     7315          solver_->setWarmStart(&bestSolutionBasis_);
    73297316      }
    73307317      // Give a hint to do dual
     
    73787365      objectiveValue = solver_->getObjValue()*solver_->getObjSense();
    73797366    }
     7367    bestSolutionBasis_ = CoinWarmStartBasis();
    73807368   
    73817369    /*
     
    75277515                        // add to global list
    75287516                        OsiRowCut newCut(*thisCut);
    7529                         newCut.setGloballyValid(2);
     7517                        newCut.setGloballyValid(true);
    75307518                        newCut.mutableRow().setTestForDuplicateIndex(false);
    75317519                        globalCuts_.insert(newCut) ;
     
    77047692              // add to global list
    77057693              OsiRowCut newCut(*thisCut);
    7706               newCut.setGloballyValid(2);
     7694              newCut.setGloballyValid(true);
    77077695              newCut.mutableRow().setTestForDuplicateIndex(false);
    77087696              globalCuts_.insert(newCut) ;
     
    77187706          // add to global list
    77197707          OsiColCut newCut(*thisCut);
    7720           newCut.setGloballyValid(2);
     7708          newCut.setGloballyValid(true);
    77217709          globalCuts_.insert(newCut) ;
    77227710        }
     
    92039191        int start = rowStart[iRow];
    92049192        thisCut.setRow(rowLength[iRow],column+start,elementByRow+start,false);
    9205         thisCut.setGloballyValid(2);
     9193        thisCut.setGloballyValid(true);
    92069194        globalCuts_.insert(thisCut) ;
    92079195      }
     
    92269214{
    92279215  OsiRowCut newCut(cut);
    9228   newCut.setGloballyValid(2);
     9216  newCut.setGloballyValid(true);
    92299217  newCut.mutableRow().setTestForDuplicateIndex(false);
    92309218  globalCuts_.insert(newCut) ;
     
    98839871  memset(bestSolution_,0,n*sizeof(double));
    98849872  memcpy(bestSolution_,solution,numberColumns*sizeof(double));
     9873}
     9874/* Do heuristics at root.
     9875   0 - don't delete
     9876   1 - delete
     9877      2 - just delete - don't even use
     9878*/
     9879void
     9880CbcModel::doHeuristicsAtRoot(int deleteHeuristicsAfterwards)
     9881{
     9882 
     9883  int numberColumns = getNumCols() ;
     9884  double * newSolution = new double [numberColumns] ;
     9885  int i;
     9886  if (deleteHeuristicsAfterwards!=2) {
     9887    if (deleteHeuristicsAfterwards) {
     9888      assert (!usedInSolution_);
     9889      usedInSolution_ = new int [numberColumns];
     9890      CoinZeroN(usedInSolution_,numberColumns);
     9891    }
     9892    double heuristicValue = getCutoff() ;
     9893    int found = -1; // no solution found
     9894    CbcEventHandler *eventHandler = getEventHandler() ;
     9895    if (eventHandler)
     9896      eventHandler->setModel(this);
     9897   
     9898    currentPassNumber_ = 1; // so root heuristics will run
     9899    for (i = 0;i<numberHeuristics_;i++) {
     9900      // see if heuristic will do anything
     9901      double saveValue = heuristicValue ;
     9902      int ifSol = heuristic_[i]->solution(heuristicValue,
     9903                                          newSolution);
     9904      if (ifSol>0) {
     9905        // better solution found
     9906        found = i ;
     9907        incrementUsed(newSolution);
     9908        // increment number of solutions so other heuristics can test
     9909        numberSolutions_++;
     9910        numberHeuristicSolutions_++;
     9911      } else {
     9912        heuristicValue = saveValue ;
     9913      }
     9914    }
     9915    currentPassNumber_ = 0;
     9916    /*
     9917      Did any of the heuristics turn up a new solution? Record it before we free
     9918      the vector.
     9919    */
     9920    if (found >= 0) {
     9921      // For compiler error on terra cluster!
     9922      if (found<numberHeuristics_)
     9923        lastHeuristic_ = heuristic_[found];
     9924      else
     9925        lastHeuristic_ = heuristic_[0];
     9926      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
     9927      CbcTreeLocal * tree
     9928        = dynamic_cast<CbcTreeLocal *> (tree_);
     9929      if (tree)
     9930        tree->passInSolution(bestSolution_,heuristicValue);
     9931      if (eventHandler) {
     9932        if (!eventHandler->event(CbcEventHandler::solution)) {
     9933          eventHappened_=true; // exit
     9934        }
     9935      }
     9936    }
     9937  }
     9938  if (!deleteHeuristicsAfterwards) {
     9939    for (i = 0;i<numberHeuristics_;i++) {
     9940      // delete FPump
     9941      CbcHeuristicFPump * pump
     9942        = dynamic_cast<CbcHeuristicFPump *> (heuristic_[i]);
     9943      if (pump) {
     9944        delete pump;
     9945        numberHeuristics_ --;
     9946        for (int j=i;j<numberHeuristics_;j++)
     9947          heuristic_[j] = heuristic_[j+1];
     9948      }
     9949    }
     9950  } else {
     9951    // delete all
     9952    for (i = 0;i<numberHeuristics_;i++)
     9953      delete heuristic_[i];
     9954    numberHeuristics_=0;
     9955    delete [] heuristic_;
     9956    heuristic_=NULL;
     9957    delete [] usedInSolution_;
     9958    usedInSolution_ = NULL;
     9959  }
     9960  delete [] newSolution ;
     9961}
     9962// Zap integer information in problem (may leave object info)
     9963void
     9964CbcModel::zapIntegerInformation(bool leaveObjects)
     9965{
     9966  numberIntegers_ = 0;
     9967  delete [] integerVariable_;
     9968  integerVariable_ = NULL;
     9969  if (!leaveObjects&&ownObjects_) {
     9970    int i;
     9971    for (i=0;i<numberObjects_;i++)
     9972      delete object_[i];
     9973    delete [] object_;
     9974    numberObjects_=0;
     9975    object_=NULL;
     9976  }
    98859977}
    98869978// Create C++ lines to get to current state
     
    1005110143  }
    1005210144}
     10145// Set original columns as created by preprocessing
     10146void
     10147CbcModel::setOriginalColumns(const int * originalColumns)
     10148{
     10149  int numberColumns = getNumCols();
     10150  delete [] originalColumns_;
     10151  originalColumns_ = CoinCopyOfArray(originalColumns,numberColumns);
     10152}
    1005310153// Set the cut modifier method
    1005410154void
  • branches/BSP/trunk/Cbc/src/CbcModel.hpp

    r742 r764  
    732732  inline int howOftenGlobalScan() const
    733733  { return howOftenGlobalScan_;}
    734   /// Original columns as created by integerPresolve
     734  /// Original columns as created by integerPresolve or preprocessing
    735735  inline int * originalColumns() const
    736736  { return originalColumns_;}
     737  /// Set original columns as created by preprocessing
     738  void setOriginalColumns(const int * originalColumns) ;
    737739
    738740  /** Set the print frequency.
     
    11041106  /// Get number of heuristic solutions
    11051107  inline int getNumberHeuristicSolutions() const { return numberHeuristicSolutions_;}
     1108  /// Set number of heuristic solutions
     1109  inline void setNumberHeuristicSolutions(int value) { numberHeuristicSolutions_=value;}
    11061110
    11071111  /// Set objective function sense (1 for min (default), -1 for max,)
     
    15051509    { return continuousSolver_;}
    15061510
     1511    /// Create solver with continuous state
     1512    inline void createContinuousSolver()
     1513    { continuousSolver_ = solver_->clone();}
     1514    /// Clear solver with continuous state
     1515    inline void clearContinuousSolver()
     1516    { delete continuousSolver_; continuousSolver_ = NULL;}
     1517
    15071518  /// A copy of the solver, taken at constructor or by saveReferenceSolver
    15081519  inline OsiSolverInterface * referenceSolver() const
     
    16171628  */
    16181629  void convertToDynamic();
     1630  /// Zap integer information in problem (may leave object info)
     1631  void zapIntegerInformation(bool leaveObjects=true);
    16191632  /// Use cliques for pseudocost information - return nonzero if infeasible
    16201633  int cliquePseudoCosts(int doStatistics);
    16211634  /// Fill in useful estimates
    16221635  void pseudoShadow(double * down, double * up);
    1623 
     1636  /** Do heuristics at root.
     1637      0 - don't delete
     1638      1 - delete
     1639      2 - just delete - don't even use
     1640  */
     1641  void doHeuristicsAtRoot(int deleteHeuristicsAfterwards=0);
    16241642  /// Get the hotstart solution
    16251643  inline const double * hotstartSolution() const
     
    16591677  /// Generate an OsiBranchingInformation object
    16601678  OsiBranchingInformation usefulInformation() const;
     1679  /** Warm start object produced by heuristic or strong branching
     1680
     1681      If get a valid integer solution outside branch and bound then it can take
     1682      a reasonable time to solve LP which produces clean solution.  If this object has
     1683      any size then it will be used in solve.
     1684  */
     1685  inline void setBestSolutionBasis(const CoinWarmStartBasis & bestSolutionBasis)
     1686  { bestSolutionBasis_ = bestSolutionBasis;}
    16611687  //@}
    16621688
     
    17341760  */
    17351761  mutable const double * testSolution_;
     1762  /** Warm start object produced by heuristic or strong branching
     1763
     1764      If get a valid integer solution outside branch and bound then it can take
     1765      a reasonable time to solve LP which produces clean solution.  If this object has
     1766      any size then it will be used in solve.
     1767  */
     1768  CoinWarmStartBasis bestSolutionBasis_ ;
    17361769  /// Global cuts
    17371770  OsiCuts globalCuts_;
     
    19511984  bool ownObjects_;
    19521985 
    1953   /// Original columns as created by integerPresolve
     1986  /// Original columns as created by integerPresolve or preprocessing
    19541987  int * originalColumns_;
    19551988  /// How often to scan global cuts
     
    20682101int callCbc1(const char * input2, CbcModel & babSolver);
    20692102int callCbc1(const std::string input2, CbcModel & babSolver);
     2103// And when CbcMain0 already called to initialize (with call back) (see CbcMain1 for whereFrom)
     2104int callCbc1(const char * input2, CbcModel & babSolver, int (CbcModel * currentSolver, int whereFrom));
     2105int callCbc1(const std::string input2, CbcModel & babSolver, int (CbcModel * currentSolver, int whereFrom));
     2106int CbcMain1 (int argc, const char *argv[],CbcModel & babSolver, int (CbcModel * currentSolver, int whereFrom), int call_CbcClpUnitTest_on_777 = 0);
    20702107#endif
  • branches/BSP/trunk/Cbc/src/CbcSolver.cpp

    r763 r764  
    3838#include "OsiRowCutDebugger.hpp"
    3939#include "OsiChooseVariable.hpp"
     40#include "OsiAuxInfo.hpp"
    4041//#define USER_HAS_FAKE_CLP
    4142//#define USER_HAS_FAKE_CBC
     
    8081}
    8182#endif
     83//#define DMALLOC
    8284#ifdef DMALLOC
    8385#include "dmalloc.h"
     
    175177   }
    176178}
     179//#define CBC_SIG_TRAP
     180#ifdef CBC_SIG_TRAP
     181#include <setjmp.h>
     182static sigjmp_buf cbc_seg_buffer;
     183extern "C" {
     184   static void signal_handler_error(int whichSignal)
     185   {
     186     siglongjmp(cbc_seg_buffer,1);
     187   }
     188}
     189#endif
    177190
    178191int CbcOrClpRead_mode=1;
     
    479492  }
    480493}
     494#if 1
     495/*
     496  On input
     497  doAction - 0 just fix in original and return NULL
     498             1 return fixed non-presolved solver
     499             2 return fixed presolved solver
     500             3 do heuristics and set best solution
     501             4 do BAB and just set best solution
     502             -2 cleanup afterwards if using 2
     503  On ouput - number fixed
     504*/
     505static OsiClpSolverInterface * fixVubs(CbcModel & model, int numberLevels,
     506                                       int & doAction, CoinMessageHandler * generalMessageHandler)
     507{
     508  OsiSolverInterface * solver = model.solver()->clone();
     509  OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
     510  ClpSimplex * lpSolver = clpSolver->getModelPtr();
     511  // Tighten bounds
     512  lpSolver->tightenPrimalBounds(0.0,11,true);
     513  char generalPrint[200];
     514  //const double *objective = solver->getObjCoefficients() ;
     515  double *columnLower = lpSolver->columnLower() ;
     516  double *columnUpper = lpSolver->columnUpper() ;
     517  int numberColumns = solver->getNumCols() ;
     518  int numberRows = solver->getNumRows();
     519  int iRow,iColumn;
     520
     521  // Row copy
     522  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
     523  const double * elementByRow = matrixByRow.getElements();
     524  const int * column = matrixByRow.getIndices();
     525  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
     526  const int * rowLength = matrixByRow.getVectorLengths();
     527
     528  // Column copy
     529  CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
     530  //const double * element = matrixByCol.getElements();
     531  const int * row = matrixByCol.getIndices();
     532  const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
     533  const int * columnLength = matrixByCol.getVectorLengths();
     534
     535  const double * rowLower = solver->getRowLower();
     536  const double * rowUpper = solver->getRowUpper();
     537
     538  // Get maximum size of VUB tree
     539  // otherColumn is one fixed to 0 if this one zero
     540  int nEl = matrixByCol.getNumElements();
     541  CoinBigIndex * fixColumn = new CoinBigIndex [numberColumns+1];
     542  int * otherColumn = new int [nEl];
     543  int * fix = new int[numberColumns];
     544  char * mark = new char [numberColumns];
     545  memset(mark,0,numberColumns);
     546  int numberInteger=0;
     547  int numberOther=0;
     548  fixColumn[0]=0;
     549  double large = lpSolver->largeValue(); // treat bounds > this as infinite
     550#ifndef NDEBUG
     551  double large2= 1.0e10*large;
     552#endif
     553  double tolerance = lpSolver->primalTolerance();
     554  int * check = new int[numberRows];
     555  for (iRow=0;iRow<numberRows;iRow++)
     556    check[iRow]=-1;
     557  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     558    fix[iColumn]=-1;
     559    if (columnUpper[iColumn] > columnLower[iColumn]+1.0e-8) {
     560      if (solver->isInteger(iColumn))
     561        numberInteger++;
     562      if (columnLower[iColumn]==0.0) {
     563        bool infeasible=false;
     564        fix[iColumn]=0;
     565        // fake upper bound
     566        double saveUpper = columnUpper[iColumn];
     567        columnUpper[iColumn]=0.0;
     568        for (CoinBigIndex i=columnStart[iColumn];
     569             i<columnStart[iColumn]+columnLength[iColumn];i++) {
     570          iRow = row[i];
     571          if (rowLower[iRow]<-1.0e6&&rowUpper[iRow]>1.0e6)
     572            continue; // unlikely
     573          //==
     574          // possible row
     575          int infiniteUpper = 0;
     576          int infiniteLower = 0;
     577          double maximumUp = 0.0;
     578          double maximumDown = 0.0;
     579          double newBound;
     580          CoinBigIndex rStart = rowStart[iRow];
     581          CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
     582          CoinBigIndex j;
     583          int kColumn;
     584          // Compute possible lower and upper ranges
     585          for (j = rStart; j < rEnd; ++j) {
     586            double value=elementByRow[j];
     587            kColumn = column[j];
     588            if (value > 0.0) {
     589              if (columnUpper[kColumn] >= large) {
     590                ++infiniteUpper;
     591              } else {
     592                maximumUp += columnUpper[kColumn] * value;
     593              }
     594              if (columnLower[kColumn] <= -large) {
     595                ++infiniteLower;
     596              } else {
     597                maximumDown += columnLower[kColumn] * value;
     598              }
     599            } else if (value<0.0) {
     600              if (columnUpper[kColumn] >= large) {
     601                ++infiniteLower;
     602              } else {
     603                maximumDown += columnUpper[kColumn] * value;
     604              }
     605              if (columnLower[kColumn] <= -large) {
     606                ++infiniteUpper;
     607              } else {
     608                maximumUp += columnLower[kColumn] * value;
     609              }
     610            }
     611          }
     612          // Build in a margin of error
     613          maximumUp += 1.0e-8*fabs(maximumUp);
     614          maximumDown -= 1.0e-8*fabs(maximumDown);
     615          double maxUp = maximumUp+infiniteUpper*1.0e31;
     616          double maxDown = maximumDown-infiniteLower*1.0e31;
     617          if (maxUp <= rowUpper[iRow] + tolerance &&
     618              maxDown >= rowLower[iRow] - tolerance) {
     619            //printf("Redundant row in vubs %d\n",iRow);
     620          } else {
     621            if (maxUp < rowLower[iRow] -100.0*tolerance ||
     622                maxDown > rowUpper[iRow]+100.0*tolerance) {
     623              infeasible=true;
     624              break;
     625            }
     626            double lower = rowLower[iRow];
     627            double upper = rowUpper[iRow];
     628            for (j = rStart; j < rEnd; ++j) {
     629              double value=elementByRow[j];
     630              kColumn = column[j];
     631              double nowLower = columnLower[kColumn];
     632              double nowUpper = columnUpper[kColumn];
     633              if (value > 0.0) {
     634                // positive value
     635                if (lower>-large) {
     636                  if (!infiniteUpper) {
     637                    assert(nowUpper < large2);
     638                    newBound = nowUpper +
     639                      (lower - maximumUp) / value;
     640                    // relax if original was large
     641                    if (fabs(maximumUp)>1.0e8)
     642                      newBound -= 1.0e-12*fabs(maximumUp);
     643                  } else if (infiniteUpper==1&&nowUpper>large) {
     644                    newBound = (lower -maximumUp) / value;
     645                    // relax if original was large
     646                    if (fabs(maximumUp)>1.0e8)
     647                      newBound -= 1.0e-12*fabs(maximumUp);
     648                  } else {
     649                    newBound = -COIN_DBL_MAX;
     650                  }
     651                  if (newBound > nowLower + 1.0e-12&&newBound>-large) {
     652                    // Tighten the lower bound
     653                    // check infeasible (relaxed)
     654                    if (nowUpper < newBound) {
     655                      if (nowUpper - newBound <
     656                          -100.0*tolerance) {
     657                        infeasible=true;
     658                        break;
     659                      }
     660                    }
     661                  }
     662                }
     663                if (upper <large) {
     664                  if (!infiniteLower) {
     665                    assert(nowLower >- large2);
     666                    newBound = nowLower +
     667                      (upper - maximumDown) / value;
     668                    // relax if original was large
     669                    if (fabs(maximumDown)>1.0e8)
     670                    newBound += 1.0e-12*fabs(maximumDown);
     671                  } else if (infiniteLower==1&&nowLower<-large) {
     672                    newBound =   (upper - maximumDown) / value;
     673                    // relax if original was large
     674                    if (fabs(maximumDown)>1.0e8)
     675                      newBound += 1.0e-12*fabs(maximumDown);
     676                  } else {
     677                    newBound = COIN_DBL_MAX;
     678                  }
     679                  if (newBound < nowUpper - 1.0e-12&&newBound<large) {
     680                    // Tighten the upper bound
     681                    // check infeasible (relaxed)
     682                    if (nowLower > newBound) {
     683                      if (newBound - nowLower <
     684                          -100.0*tolerance) {
     685                        infeasible=true;
     686                        break;
     687                      } else {
     688                        newBound=nowLower;
     689                      }
     690                    }
     691                    if (!newBound||(solver->isInteger(kColumn)&&newBound<0.999)) {
     692                      // fix to zero
     693                      if (!mark[kColumn]) {
     694                        otherColumn[numberOther++]=kColumn;
     695                        mark[kColumn]=1;
     696                        if (check[iRow]==-1)
     697                          check[iRow]=iColumn;
     698                        else
     699                          assert(check[iRow]==iColumn);
     700                      }
     701                    }
     702                  }
     703                }
     704              } else {
     705                // negative value
     706                if (lower>-large) {
     707                  if (!infiniteUpper) {
     708                    assert(nowLower < large2);
     709                    newBound = nowLower +
     710                      (lower - maximumUp) / value;
     711                    // relax if original was large
     712                    if (fabs(maximumUp)>1.0e8)
     713                      newBound += 1.0e-12*fabs(maximumUp);
     714                  } else if (infiniteUpper==1&&nowLower<-large) {
     715                    newBound = (lower -maximumUp) / value;
     716                    // relax if original was large
     717                    if (fabs(maximumUp)>1.0e8)
     718                      newBound += 1.0e-12*fabs(maximumUp);
     719                  } else {
     720                    newBound = COIN_DBL_MAX;
     721                  }
     722                  if (newBound < nowUpper - 1.0e-12&&newBound<large) {
     723                    // Tighten the upper bound
     724                    // check infeasible (relaxed)
     725                    if (nowLower > newBound) {
     726                      if (newBound - nowLower <
     727                          -100.0*tolerance) {
     728                        infeasible=true;
     729                        break;
     730                      } else {
     731                        newBound=nowLower;
     732                      }
     733                    }
     734                    if (!newBound||(solver->isInteger(kColumn)&&newBound<0.999)) {
     735                      // fix to zero
     736                      if (!mark[kColumn]) {
     737                        otherColumn[numberOther++]=kColumn;
     738                        mark[kColumn]=1;
     739                        if (check[iRow]==-1)
     740                          check[iRow]=iColumn;
     741                        else
     742                          assert(check[iRow]==iColumn);
     743                      }
     744                    }
     745                  }
     746                }
     747                if (upper <large) {
     748                  if (!infiniteLower) {
     749                    assert(nowUpper < large2);
     750                    newBound = nowUpper +
     751                      (upper - maximumDown) / value;
     752                    // relax if original was large
     753                    if (fabs(maximumDown)>1.0e8)
     754                      newBound -= 1.0e-12*fabs(maximumDown);
     755                  } else if (infiniteLower==1&&nowUpper>large) {
     756                    newBound =   (upper - maximumDown) / value;
     757                    // relax if original was large
     758                    if (fabs(maximumDown)>1.0e8)
     759                      newBound -= 1.0e-12*fabs(maximumDown);
     760                  } else {
     761                    newBound = -COIN_DBL_MAX;
     762                  }
     763                  if (newBound > nowLower + 1.0e-12&&newBound>-large) {
     764                    // Tighten the lower bound
     765                    // check infeasible (relaxed)
     766                    if (nowUpper < newBound) {
     767                      if (nowUpper - newBound <
     768                          -100.0*tolerance) {
     769                        infeasible=true;
     770                        break;
     771                      }
     772                    }
     773                  }
     774                }
     775              }
     776            }
     777          }
     778        }
     779        for (int i=fixColumn[iColumn];i<numberOther;i++)
     780          mark[otherColumn[i]]=0;
     781        // reset bound unless infeasible
     782        if (!infeasible||!solver->isInteger(iColumn))
     783          columnUpper[iColumn]=saveUpper;
     784        else if (solver->isInteger(iColumn))
     785          columnLower[iColumn]=1.0;
     786      }
     787    }
     788    fixColumn[iColumn+1]=numberOther;
     789  }
     790  delete [] check;
     791  delete [] mark;
     792  // Now do reverse way
     793  int * counts = new int [numberColumns];
     794  CoinZeroN(counts,numberColumns);
     795  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     796    for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++)
     797      counts[otherColumn[i]]++;
     798  }
     799  numberOther=0;
     800  CoinBigIndex * fixColumn2 = new CoinBigIndex [numberColumns+1];
     801  int * otherColumn2 = new int [fixColumn[numberColumns]];
     802  fixColumn2[0]=0;
     803  for ( iColumn=0;iColumn<numberColumns;iColumn++) {
     804    numberOther += counts[iColumn];
     805    counts[iColumn]=0;
     806    fixColumn2[iColumn+1]=numberOther;
     807  }
     808  for ( iColumn=0;iColumn<numberColumns;iColumn++) {
     809    for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
     810      int jColumn=otherColumn[i];
     811      CoinBigIndex put = fixColumn2[jColumn] + counts[jColumn];
     812      counts[jColumn]++;
     813      otherColumn2[put]=iColumn;
     814    }
     815  }
     816  // get top layer i.e. those which are not fixed by any other
     817  int kLayer=0;
     818  while (true) {
     819    int numberLayered=0;
     820    int numberInteger=0;
     821    for ( iColumn=0;iColumn<numberColumns;iColumn++) {
     822      if (fix[iColumn]==kLayer) {
     823        for (int i=fixColumn2[iColumn];i<fixColumn2[iColumn+1];i++) {
     824          int jColumn=otherColumn2[i];
     825          if (fix[jColumn]==kLayer) {
     826            fix[iColumn]=kLayer+1;
     827          }
     828        }
     829      }
     830      if (fix[iColumn]==kLayer) {
     831        numberLayered++;
     832        if (solver->isInteger(iColumn))
     833          numberInteger++;
     834      }
     835    }
     836    if (numberLayered) {
     837      printf("%d (%d integer) at priority %d\n",numberLayered,numberInteger,kLayer);
     838      kLayer++;
     839    } else {
     840      break;
     841    }
     842  }
     843  delete [] counts;
     844  delete [] fixColumn;
     845  delete [] otherColumn;
     846  delete [] otherColumn2;
     847  delete [] fixColumn2;
     848  // Now do fixing
     849  delete [] fix;
     850  delete solver;
     851  return NULL;
     852}
     853#endif
    481854#ifdef COIN_HAS_ASL
    482855/*  Returns OsiSolverInterface (User should delete)
     
    11851558#endif
    11861559}
    1187 int callCbc1(const char * input2, CbcModel & model)
     1560int callCbc1(const char * input2, CbcModel & model, int callBack(CbcModel * currentSolver, int whereFrom))
    11881561{
    11891562  char * input = strdup(input2);
     
    12361609  CbcOrClpReadCommand=stdin;
    12371610  noPrinting=false;
    1238   int returnCode = CbcMain1(n+2,const_cast<const char **>(argv),model);
     1611  int returnCode = CbcMain1(n+2,const_cast<const char **>(argv),model,callBack);
    12391612  for (int k=0;k<n+2;k++)
    12401613    free(argv[k]);
     
    12891662  return returnCode;
    12901663}
    1291 
     1664static int dummyCallBack(CbcModel * model, int whereFrom)
     1665{
     1666  return 0;
     1667}
     1668int CbcMain1 (int argc, const char *argv[],
     1669              CbcModel  & model, int call_CbcClpUnitTest_on_777)
     1670{
     1671  return CbcMain1(argc,argv,model,dummyCallBack, call_CbcClpUnitTest_on_777);
     1672}
     1673int callCbc1(const std::string input2, CbcModel & babSolver, int callBack(CbcModel * currentSolver, int whereFrom))
     1674{
     1675  char * input3 = strdup(input2.c_str());
     1676  int returnCode=callCbc1(input3,babSolver,callBack);
     1677  free(input3);
     1678  return returnCode;
     1679}
     1680int callCbc1(const char * input2, CbcModel & model)
     1681{
     1682  return callCbc1(input2,model,dummyCallBack);
     1683}
    12921684int CbcMain (int argc, const char *argv[],
    12931685             CbcModel  & model)
     
    14221814  parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption("on");
    14231815  parameters[whichParam(RINS,numberParameters,parameters)].setCurrentOption("off");
     1816  parameters[whichParam(LOCALTREE,numberParameters,parameters)].setCurrentOption("off");
    14241817  parameters[whichParam(COSTSTRATEGY,numberParameters,parameters)].setCurrentOption("off");
    14251818}
     1819
     1820/* 1 - add heuristics to model
     1821   2 - do heuristics (and set cutoff and best solution)
     1822   3 - for miplib test so skip some
     1823*/
     1824static int doHeuristics(CbcModel * model,int type)
     1825{
     1826  bool anyToDo=false;
     1827  int logLevel = parameters[whichParam(LOGLEVEL,numberParameters,parameters)].intValue();
     1828  int useFpump = parameters[whichParam(FPUMP,numberParameters,parameters)].currentOptionAsInteger();
     1829  int useRounding = parameters[whichParam(ROUNDING,numberParameters,parameters)].currentOptionAsInteger();
     1830  int useGreedy = parameters[whichParam(GREEDY,numberParameters,parameters)].currentOptionAsInteger();
     1831  int useCombine = parameters[whichParam(COMBINE,numberParameters,parameters)].currentOptionAsInteger();
     1832  int useRINS = parameters[whichParam(RINS,numberParameters,parameters)].currentOptionAsInteger();
     1833  // FPump done first as it only works if no solution
     1834  int kType = (type<3) ? type : 1;
     1835  if (useFpump>=kType) {
     1836    anyToDo=true;
     1837    CbcHeuristicFPump heuristic4(*model);
     1838    heuristic4.setFractionSmall(0.5);
     1839    double dextra3 = parameters[whichParam(DEXTRA3,numberParameters,parameters)].doubleValue();
     1840    if (dextra3)
     1841      heuristic4.setFractionSmall(dextra3);
     1842    heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
     1843    int pumpTune=parameters[whichParam(FPUMPTUNE,numberParameters,parameters)].intValue();
     1844    if (pumpTune>0) {
     1845      /*
     1846        >=10000000 for using obj
     1847        >=1000000 use as accumulate switch
     1848        >=1000 use index+1 as number of large loops
     1849        >=100 use 0.05 objvalue as increment
     1850        >=10 use +0.1 objvalue for cutoff (add)
     1851        1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
     1852        4 and static continuous, 5 as 3 but no internal integers
     1853        6 as 3 but all slack basis!
     1854      */
     1855      double value = model->solver()->getObjSense()*model->solver()->getObjValue();
     1856      int w = pumpTune/10;
     1857      int c = w % 10;
     1858      w /= 10;
     1859      int i = w % 10;
     1860      w /= 10;
     1861      int r = w;
     1862      int accumulate = r/1000;
     1863      r -= 1000*accumulate;
     1864      if (accumulate>=10) {
     1865        int which = accumulate/10;
     1866        accumulate -= 10*which;
     1867        which--;
     1868        // weights and factors
     1869        double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
     1870        double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
     1871        heuristic4.setInitialWeight(weight[which]);
     1872        heuristic4.setWeightFactor(factor[which]);
     1873      }
     1874      // fake cutoff
     1875      if (logLevel>1)
     1876        printf("Setting ");
     1877      if (c) {
     1878        double cutoff;
     1879        model->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
     1880        cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
     1881        double dextra1 = parameters[whichParam(DEXTRA1,numberParameters,parameters)].doubleValue();
     1882        if (dextra1)
     1883          cutoff=dextra1;
     1884        heuristic4.setFakeCutoff(cutoff);
     1885        if (logLevel>1)
     1886          printf("fake cutoff of %g ",cutoff);
     1887      }
     1888      if (i||r) {
     1889        // also set increment
     1890        //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
     1891        double increment = 0.0;
     1892        double dextra2 = parameters[whichParam(DEXTRA2,numberParameters,parameters)].doubleValue();
     1893        if (dextra2)
     1894          increment = dextra2;
     1895        heuristic4.setAbsoluteIncrement(increment);
     1896        heuristic4.setAccumulate(accumulate);
     1897        heuristic4.setMaximumRetries(r+1);
     1898        if (logLevel>1) {
     1899          if (i)
     1900            printf("increment of %g ",heuristic4.absoluteIncrement());
     1901          if (accumulate)
     1902            printf("accumulate of %d ",accumulate);
     1903          printf("%d retries ",r+2);
     1904        }
     1905      }
     1906      pumpTune = pumpTune%100;
     1907      if (logLevel>1)
     1908        printf("and setting when to %d\n",pumpTune+10);
     1909      if (pumpTune==6)
     1910        pumpTune =13;
     1911      heuristic4.setWhen(pumpTune+10);
     1912    }
     1913    heuristic4.setHeuristicName("feasibility pump");
     1914    //#define ROLF
     1915#ifdef ROLF   
     1916    CbcHeuristicFPump pump(*model);
     1917    pump.setMaximumTime(60);
     1918    pump.setMaximumPasses(100);
     1919    pump.setMaximumRetries(1);
     1920    pump.setFixOnReducedCosts(0);
     1921    pump.setHeuristicName("Feasibility pump");
     1922    pump.setFractionSmall(1.0);
     1923    pump.setWhen(13);
     1924    model->addHeuristic(&pump);
     1925#else
     1926    model->addHeuristic(&heuristic4);
     1927#endif
     1928  }
     1929  if (useRounding>=type) {
     1930    CbcRounding heuristic1(*model);
     1931    heuristic1.setHeuristicName("rounding");
     1932    model->addHeuristic(&heuristic1) ;
     1933    anyToDo=true;
     1934  }
     1935  if (useCombine>=type) {
     1936    CbcHeuristicLocal heuristic2(*model);
     1937    heuristic2.setHeuristicName("combine solutions");
     1938    heuristic2.setFractionSmall(0.6);
     1939    heuristic2.setSearchType(1);
     1940    model->addHeuristic(&heuristic2);
     1941    anyToDo=true;
     1942  }
     1943  if (useGreedy>=type) {
     1944    CbcHeuristicGreedyCover heuristic3(*model);
     1945    heuristic3.setHeuristicName("greedy cover");
     1946    CbcHeuristicGreedyEquality heuristic3a(*model);
     1947    heuristic3a.setHeuristicName("greedy equality");
     1948    model->addHeuristic(&heuristic3);
     1949    model->addHeuristic(&heuristic3a);
     1950    anyToDo=true;
     1951  }
     1952  if (useRINS>=kType) {
     1953    CbcHeuristicRINS heuristic5(*model);
     1954    heuristic5.setHeuristicName("RINS");
     1955    heuristic5.setFractionSmall(0.6);
     1956    model->addHeuristic(&heuristic5) ;
     1957    anyToDo=true;
     1958  }
     1959  if (type==2&&anyToDo) {
     1960    // Do heuristics
     1961#if 1
     1962    // clean copy
     1963    CbcModel model2(*model);
     1964    // But get rid of heuristics in model
     1965    model->doHeuristicsAtRoot(2);
     1966    if (logLevel<=1)
     1967      model2.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
     1968    OsiBabSolver defaultC;
     1969    //solver_->setAuxiliaryInfo(&defaultC);
     1970    model2.passInSolverCharacteristics(&defaultC);
     1971    // Save bounds
     1972    int numberColumns = model2.solver()->getNumCols();
     1973    model2.createContinuousSolver();
     1974    bool cleanModel = !model2.numberIntegers()&&!model2.numberObjects();
     1975    model2.findIntegers(false);
     1976    model2.doHeuristicsAtRoot(1);
     1977    if (cleanModel)
     1978      model2.zapIntegerInformation(false);
     1979    if (model2.bestSolution()) {
     1980      double value = model2.getMinimizationObjValue();
     1981      model->setCutoff(value);
     1982      model->setBestSolution(model2.bestSolution(),numberColumns,value);
     1983      model->setSolutionCount(1);
     1984      model->setNumberHeuristicSolutions(1);
     1985    }
     1986#else
     1987    if (logLevel<=1)
     1988      model->solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
     1989    OsiBabSolver defaultC;
     1990    //solver_->setAuxiliaryInfo(&defaultC);
     1991    model->passInSolverCharacteristics(&defaultC);
     1992    // Save bounds
     1993    int numberColumns = model->solver()->getNumCols();
     1994    model->createContinuousSolver();
     1995    bool cleanModel = !model->numberIntegers()&&!model->numberObjects();
     1996    model->findIntegers(false);
     1997    model->doHeuristicsAtRoot(true);
     1998    if (cleanModel)
     1999      model->zapIntegerInformation(false);
     2000#endif
     2001    return 0;
     2002  } else {
     2003    return 0;
     2004  }
     2005}
     2006
     2007void CbcClpUnitTest (const CbcModel & saveModel,
     2008                     std::string& dirMiplib, bool unitTestOnly);
     2009
     2010/* Meaning of whereFrom:
     2011   1 after initial solve by dualsimplex etc
     2012   2 after preprocessing
     2013   3 just before branchAndBound (so user can override)
     2014   4 just after branchAndBound (before postprocessing)
     2015   5 after postprocessing
     2016   6 after a user called heuristic phase
     2017*/
     2018
    14262019int CbcMain1 (int argc, const char *argv[],
    1427               CbcModel  & model, int call_CbcClpUnitTest_on_777)
     2020              CbcModel  & model,
     2021              int callBack(CbcModel * currentSolver, int whereFrom),
     2022              int call_CbcClpUnitTest_on_777)
    14282023{
    14292024  /* Note
     
    14602055  }
    14612056  CbcModel * babModel = NULL;
     2057  double time0;
    14622058  {
    14632059    double time1 = CoinCpuTime(),time2;
     2060    time0=time1;
    14642061    bool goodModel=(originalSolver->getNumCols()) ? true : false;
    14652062
     
    18092406    bool storedCuts = false;
    18102407
    1811     bool useRounding=true;
    1812     bool useFpump=true;
    1813     bool useGreedy=true;
    1814     bool useCombine=true;
    1815     bool useLocalTree=false;
    1816     bool useRINS=false;
    18172408    int useCosts=0;
    18182409    // don't use input solution
     
    18552446      // next command
    18562447      field=CoinReadGetCommand(argc,argv);
     2448      // Reset time
     2449      time1 = CoinCpuTime();
    18572450      // adjust field if has odd trailing characters
    18582451      char temp [200];
     
    19272520          }
    19282521#endif
     2522          lpSolver = clpSolver->getModelPtr();
    19292523          if (!lpSolver->integerInformation()&&!numberSOS&&
    19302524              !clpSolver->numberSOS()&&!model.numberObjects()&&!clpSolver->numberObjects())
     
    20612655              case DJFIX:
    20622656                djFix=value;
    2063                 preSolve=5;
    2064                 defaultSettings=false; // user knows what she is doing
     2657                if (goodModel&&djFix<1.0e20) {
     2658                  // do some fixing
     2659                  clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
     2660                  clpSolver->initialSolve();
     2661                  lpSolver = clpSolver->getModelPtr();
     2662                  int numberColumns = lpSolver->numberColumns();
     2663                  int i;
     2664                  const char * type = lpSolver->integerInformation();
     2665                  double * lower = lpSolver->columnLower();
     2666                  double * upper = lpSolver->columnUpper();
     2667                  double * solution = lpSolver->primalColumnSolution();
     2668                  double * dj = lpSolver->dualColumnSolution();
     2669                  int numberFixed=0;
     2670                  double dextra4 = parameters[whichParam(DEXTRA4,numberParameters,parameters)].doubleValue();
     2671                  if (dextra4)
     2672                    printf("Multiple for continuous dj fixing is %g\n",dextra4);
     2673                  for (i=0;i<numberColumns;i++) {
     2674                    double djValue = dj[i];
     2675                    if (!type[i])
     2676                      djValue *= dextra4;
     2677                    if (type[i]||dextra4) {
     2678                      double value = solution[i];
     2679                      if (value<lower[i]+1.0e-5&&djValue>djFix) {
     2680                        solution[i]=lower[i];
     2681                        upper[i]=lower[i];
     2682                        numberFixed++;
     2683                      } else if (value>upper[i]-1.0e-5&&djValue<-djFix) {
     2684                        solution[i]=upper[i];
     2685                        lower[i]=upper[i];
     2686                        numberFixed++;
     2687                      }
     2688                    }
     2689                  }
     2690                  sprintf(generalPrint,"%d columns fixed\n",numberFixed);
     2691                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
     2692                    << generalPrint
     2693                    <<CoinMessageEol;
     2694                }
    20652695                break;
    20662696              case GAPRATIO:
     
    20772707            }
    20782708          } else if (valid==1) {
    2079             abort();
     2709            std::cout<<" is illegal for double parameter "<<parameters[iParam].name()<<" value remains "<<
     2710              parameters[iParam].doubleValue()<<std::endl;
    20802711          } else {
    20812712            std::cout<<parameters[iParam].name()<<" has value "<<
     
    21172748              else if (parameters[iParam].type()==CUTPASSINTREE)
    21182749                cutPassInTree = value;
    2119               else if (parameters[iParam].type()==FPUMPITS)
    2120                 { useFpump = true;parameters[iParam].setIntValue(value);}
    21212750              else if (parameters[iParam].type()==STRONGBRANCHING||
    21222751                       parameters[iParam].type()==NUMBERBEFORE)
     
    21252754            }
    21262755          } else if (valid==1) {
    2127             abort();
     2756            std::cout<<" is illegal for integer parameter "<<parameters[iParam].name()<<" value remains "<<
     2757              parameters[iParam].intValue()<<std::endl;
    21282758          } else {
    21292759            std::cout<<parameters[iParam].name()<<" has value "<<
     
    23122942            case ROUNDING:
    23132943              defaultSettings=false; // user knows what she is doing
    2314               useRounding = (action!=0);
    23152944              break;
    23162945            case FPUMP:
    23172946              defaultSettings=false; // user knows what she is doing
    2318               useFpump=(action!=0);
    23192947              break;
    23202948            case RINS:
    2321               useRINS=(action!=0);
    23222949              break;
    23232950            case CUTSSTRATEGY:
     
    23472974              break;
    23482975            case HEURISTICSTRATEGY:
    2349               useRounding = (action!=0);
    2350               useGreedy = (action!=0);
    2351               useCombine = (action!=0);
    2352               //useLocalTree = action;
    2353               useFpump=(action!=0);
    23542976              parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption(action);
    23552977              parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption(action);
     
    23602982            case GREEDY:
    23612983              defaultSettings=false; // user knows what she is doing
    2362               useGreedy = (action!=0);
    23632984              break;
    23642985            case COMBINE:
    23652986              defaultSettings=false; // user knows what she is doing
    2366               useCombine = (action!=0);
    23672987              break;
    23682988            case LOCALTREE:
    23692989              defaultSettings=false; // user knows what she is doing
    2370               useLocalTree = (action!=0);
    23712990              break;
    23722991            case COSTSTRATEGY:
     
    24383057                //printf("dualize %d\n",dualize);
    24393058                model2 = ((ClpSimplexOther *) model2)->dualOfModel();
    2440                 sprintf(generalPrint,"Dual of model has %d rows and %d columns\n",
     3059                sprintf(generalPrint,"Dual of model has %d rows and %d columns",
    24413060                       model2->numberRows(),model2->numberColumns());
    24423061                generalMessageHandler->message(CLP_GENERAL,generalMessages)
     
    24563075                  if (preSolve<=100) {
    24573076                    presolveType=ClpSolve::presolveNumber;
    2458                     sprintf(generalPrint,"Doing %d presolve passes - picking up non-costed slacks\n",
     3077                    sprintf(generalPrint,"Doing %d presolve passes - picking up non-costed slacks",
    24593078                           preSolve);
    24603079                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
     
    24653084                    preSolve -=100;
    24663085                    presolveType=ClpSolve::presolveNumberCost;
    2467                     sprintf(generalPrint,"Doing %d presolve passes - picking up costed slacks\n",
     3086                    sprintf(generalPrint,"Doing %d presolve passes - picking up costed slacks",
    24683087                           preSolve);
    24693088                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
     
    26583277                  babModel->setProblemStatus(iStatus);
    26593278                  babModel->setSecondaryStatus(iStatus2);
     3279                }
     3280                int returnCode=callBack(&model,1);
     3281                if (returnCode) {
     3282                  // exit if user wants
     3283                  delete babModel;
     3284                  return returnCode;
    26603285                }
    26613286              }
     
    28073432              int nOut = numberRows-clpSolver->getNumRows();
    28083433              if (nOut&&!noPrinting)
    2809                 sprintf(generalPrint,"%d rows eliminated\n",nOut);
     3434                sprintf(generalPrint,"%d rows eliminated",nOut);
    28103435              generalMessageHandler->message(CLP_GENERAL,generalMessages)
    28113436                << generalPrint
     
    28343459            } else {
    28353460              std::cout<<"** Current model not valid"<<std::endl;
     3461            }
     3462            break;
     3463          case DOHEURISTIC:
     3464            if (goodModel) {
     3465              int vubAction = parameters[whichParam(VUBTRY,numberParameters,parameters)].intValue();
     3466              if (vubAction!=-1) {
     3467                // look at vubs
     3468                OsiClpSolverInterface * newSolver = fixVubs(model,1,vubAction,generalMessageHandler);
     3469                assert (!newSolver);
     3470              }
     3471              // Actually do heuristics
     3472              doHeuristics(&model,2);
     3473              if (model.bestSolution()) {
     3474                model.setProblemStatus(1);
     3475                model.setSecondaryStatus(6);
     3476              }
     3477              int returnCode=callBack(&model,6);
     3478              if (returnCode) {
     3479                // exit if user wants
     3480                delete babModel;
     3481                return returnCode;
     3482              }
    28363483            }
    28373484            break;
     
    30743721                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
    30753722                ClpSimplex * clpSolver = si->getModelPtr();
     3723                int iStatus = clpSolver->status();
     3724                int iStatus2 = clpSolver->secondaryStatus();
     3725                if (iStatus==0) {
     3726                  iStatus2=0;
     3727                } else if (iStatus==1) {
     3728                  iStatus=0;
     3729                  iStatus2=1; // say infeasible
     3730                } else if (iStatus==2) {
     3731                  iStatus=0;
     3732                  iStatus2=7; // say unbounded
     3733                } else if (iStatus==3) {
     3734                  iStatus=1;
     3735                  if (iStatus2==9)
     3736                    iStatus2=4;
     3737                  else
     3738                    iStatus2=3; // Use nodes - as closer than solutions
     3739                } else if (iStatus==4) {
     3740                  iStatus=2; // difficulties
     3741                  iStatus2=0;
     3742                }
     3743                model.setProblemStatus(iStatus);
     3744                model.setSecondaryStatus(iStatus2);
     3745                si->setWarmStart(NULL);
     3746                int returnCode=callBack(&model,1);
     3747                if (returnCode) {
     3748                  // exit if user wants
     3749                  delete babModel;
     3750                  return returnCode;
     3751                }
    30763752                clpSolver->setSpecialOptions(clpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
    30773753                if (!noPrinting) {
     
    31583834                  }
    31593835                }
    3160                 if (djFix<1.0e20) {
    3161                   // do some fixing
    3162                   int numberColumns = modelC->numberColumns();
    3163                   int i;
    3164                   const char * type = modelC->integerInformation();
    3165                   double * lower = modelC->columnLower();
    3166                   double * upper = modelC->columnUpper();
    3167                   double * solution = modelC->primalColumnSolution();
    3168                   double * dj = modelC->dualColumnSolution();
    3169                   int numberFixed=0;
    3170                   for (i=0;i<numberColumns;i++) {
    3171                     if (type[i]) {
    3172                       double value = solution[i];
    3173                       if (value<lower[i]+1.0e-5&&dj[i]>djFix) {
    3174                         solution[i]=lower[i];
    3175                         upper[i]=lower[i];
    3176                         numberFixed++;
    3177                       } else if (value>upper[i]-1.0e-5&&dj[i]<-djFix) {
    3178                         solution[i]=upper[i];
    3179                         lower[i]=upper[i];
    3180                         numberFixed++;
    3181                       }
    3182                     }
    3183                   }
    3184                   sprintf(generalPrint,"%d columns fixed\n",numberFixed);
    3185                   generalMessageHandler->message(CLP_GENERAL,generalMessages)
    3186                     << generalPrint
    3187                     <<CoinMessageEol;
    3188                 }
    3189               }
     3836              }
    31903837              // See if we want preprocessing
    31913838              OsiSolverInterface * saveSolver=NULL;
     
    33854032                if (!noPrinting) {
    33864033                  if (!solver2) {
    3387                     sprintf(generalPrint,"Pre-processing says infeasible or unbounded\n");
     4034                    sprintf(generalPrint,"Pre-processing says infeasible or unbounded");
    33884035                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
    33894036                      << generalPrint
    33904037                      <<CoinMessageEol;
    3391                     break;
    33924038                  } else {
    33934039                    //printf("processed model has %d rows, %d columns and %d elements\n",
    33944040                    //     solver2->getNumRows(),solver2->getNumCols(),solver2->getNumElements());
    33954041                  }
    3396                 }
     4042                }
     4043                if (!solver2) {
     4044                  model.setProblemStatus(0);
     4045                  model.setSecondaryStatus(1);
     4046                  babModel->setProblemStatus(0);
     4047                  babModel->setSecondaryStatus(1);
     4048                } else {
     4049                  model.setProblemStatus(-1);
     4050                  babModel->setProblemStatus(-1);
     4051                }
     4052                int returnCode=callBack(babModel,2);
     4053                if (returnCode) {
     4054                  // exit if user wants
     4055                  delete babModel;
     4056                  return returnCode;
     4057                }
     4058                if (!solver2)
     4059                  break;
     4060                if (model.bestSolution()) {
     4061                  // need to redo - in case no better found in BAB
     4062                  // just get integer part right
     4063                  const int * originalColumns = process.originalColumns();
     4064                  int numberColumns = solver2->getNumCols();
     4065                  double * bestSolution = babModel->bestSolution();
     4066                  const double * oldBestSolution = model.bestSolution();
     4067                  for (int i=0;i<numberColumns;i++) {
     4068                    int jColumn = originalColumns[i];
     4069                    bestSolution[i]=oldBestSolution[jColumn];
     4070                  }
     4071                }
    33974072                //solver2->resolve();
    33984073                if (preProcess==2) {
     
    34054080                solver2 = solver2->clone();
    34064081                babModel->assignSolver(solver2);
     4082                babModel->setOriginalColumns(process.originalColumns());
    34074083                babModel->initialSolve();
    34084084                babModel->setMaximumSeconds(timeLeft-(CoinCpuTime()-time1));
     
    35444220                delete [] dsort;
    35454221              }
    3546               // FPump done first as it only works if no solution
    3547               CbcHeuristicFPump heuristic4(*babModel);
    3548               heuristic4.setFractionSmall(0.6);
    3549               if (useFpump) {
    3550                 heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
    3551                 int pumpTune=parameters[whichParam(FPUMPTUNE,numberParameters,parameters)].intValue();
    3552                 if (pumpTune>0) {
    3553                   /*
    3554                     >=10000000 for using obj
    3555                     >=1000000 use as accumulate switch
    3556                     >=1000 use index+1 as number of large loops
    3557                     >=100 use 0.05 objvalue as increment
    3558                     >=10 use +0.1 objvalue for cutoff (add)
    3559                     1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
    3560                     4 and static continuous, 5 as 3 but no internal integers
    3561                     6 as 3 but all slack basis!
    3562                   */
    3563                   int logLevel = parameters[log].intValue();
    3564                   double value = babModel->solver()->getObjSense()*babModel->solver()->getObjValue();
    3565                   int w = pumpTune/10;
    3566                   int c = w % 10;
    3567                   w /= 10;
    3568                   int i = w % 10;
    3569                   w /= 10;
    3570                   int r = w;
    3571                   int accumulate = r/1000;
    3572                   r -= 1000*accumulate;
    3573                   if (accumulate>=10) {
    3574                     int which = accumulate/10;
    3575                     accumulate -= 10*which;
    3576                     which--;
    3577                     // weights and factors
    3578                     double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
    3579                     double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
    3580                     heuristic4.setInitialWeight(weight[which]);
    3581                     heuristic4.setWeightFactor(factor[which]);
    3582                   }
    3583                   // fake cutoff
    3584                   if (logLevel>1)
    3585                     printf("Setting ");
    3586                   if (c) {
    3587                     double cutoff;
    3588                     babModel->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
    3589                     cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
    3590                     double dextra1 = parameters[whichParam(DEXTRA1,numberParameters,parameters)].doubleValue();
    3591                     if (dextra1)
    3592                       cutoff=dextra1;
    3593                     heuristic4.setFakeCutoff(cutoff);
    3594                     if (logLevel>1)
    3595                       printf("fake cutoff of %g ",cutoff);
    3596                   }
    3597                   if (i||r) {
    3598                     // also set increment
    3599                     double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
    3600                     double dextra2 = parameters[whichParam(DEXTRA2,numberParameters,parameters)].doubleValue();
    3601                     if (dextra2)
    3602                       increment = dextra2;
    3603                     heuristic4.setAbsoluteIncrement(increment);
    3604                     heuristic4.setAccumulate(accumulate);
    3605                     heuristic4.setMaximumRetries(r+1);
    3606                     if (logLevel>1) {
    3607                       if (i)
    3608                         printf("increment of %g ",heuristic4.absoluteIncrement());
    3609                       if (accumulate)
    3610                         printf("accumulate of %d ",accumulate);
    3611                       printf("%d retries ",r+2);
    3612                     }
    3613                   }
    3614                   pumpTune = pumpTune%100;
    3615                   if (logLevel>1)
    3616                     printf("and setting when to %d\n",pumpTune+10);
    3617                   if (pumpTune==6)
    3618                     pumpTune =13;
    3619                   heuristic4.setWhen(pumpTune+10);
    3620                 }
    3621                 heuristic4.setHeuristicName("feasibility pump");
    3622                 babModel->addHeuristic(&heuristic4);
    3623               }
     4222              // Set up heuristics
     4223              doHeuristics(babModel,(!miplib) ? 1 : 3);
    36244224              if (!miplib) {
    3625                 CbcRounding heuristic1(*babModel);
    3626                 heuristic1.setHeuristicName("rounding");
    3627                 if (useRounding)
    3628                   babModel->addHeuristic(&heuristic1) ;
    3629                 CbcHeuristicLocal heuristic2(*babModel);
    3630                 heuristic2.setHeuristicName("combine solutions");
    3631                 heuristic2.setFractionSmall(0.6);
    3632                 heuristic2.setSearchType(1);
    3633                 if (useCombine)
    3634                   babModel->addHeuristic(&heuristic2);
    3635                 CbcHeuristicGreedyCover heuristic3(*babModel);
    3636                 heuristic3.setHeuristicName("greedy cover");
    3637                 CbcHeuristicGreedyEquality heuristic3a(*babModel);
    3638                 heuristic3a.setHeuristicName("greedy equality");
    3639                 if (useGreedy) {
    3640                   babModel->addHeuristic(&heuristic3);
    3641                   babModel->addHeuristic(&heuristic3a);
    3642                 }
    3643                 if (useLocalTree) {
     4225                if(parameters[whichParam(LOCALTREE,numberParameters,parameters)].currentOptionAsInteger()) {
    36444226                  CbcTreeLocal localTree(babModel,NULL,10,0,0,10000,2000);
    36454227                  babModel->passInTreeHandler(localTree);
    36464228                }
    36474229              }
    3648               CbcHeuristicRINS heuristic5(*babModel);
    3649               heuristic5.setHeuristicName("RINS");
    3650               heuristic5.setFractionSmall(0.6);
    3651               if (useRINS)
    3652                 babModel->addHeuristic(&heuristic5) ;
    36534230              if (type==MIPLIB) {
    36544231                if (babModel->numberStrong()==5&&babModel->numberBeforeTrust()==5)
     
    37894366              int mipOptions = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()%10000;
    37904367              if (mipOptions!=(1)) {
    3791                 sprintf(generalPrint,"mip options %d\n",mipOptions);
     4368                sprintf(generalPrint,"mip options %d",mipOptions);
    37924369                generalMessageHandler->message(CLP_GENERAL,generalMessages)
    37934370                  << generalPrint
     
    38104387                int moreMipOptions = parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].intValue();
    38114388                if (moreMipOptions>=0) {
    3812                   sprintf(generalPrint,"more mip options %d\n",moreMipOptions);
     4389                  sprintf(generalPrint,"more mip options %d",moreMipOptions);
    38134390                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
    38144391                    << generalPrint
     
    39414518                      }
    39424519                      babModel->setNumberObjects(n);
     4520                      babModel->zapIntegerInformation();
    39434521                    }
    39444522                    int nMissing=0;
     
    39944572                    }
    39954573                    if (nMissing) {
    3996                       sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
     4574                      sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?",nMissing);
    39974575                      generalMessageHandler->message(CLP_GENERAL,generalMessages)
    39984576                        << generalPrint
     
    40554633                        delete [] back;
    40564634                        if (nMissing) {
    4057                           sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
     4635                          sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?",nMissing);
    40584636                          generalMessageHandler->message(CLP_GENERAL,generalMessages)
    40594637                            << generalPrint
     
    42244802                        delete [] back;
    42254803                        if (nMissing) {
    4226                           sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
     4804                          sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?",nMissing);
    42274805                          generalMessageHandler->message(CLP_GENERAL,generalMessages)
    42284806                            << generalPrint
     
    43964974                }
    43974975                if (testOsiOptions>=0) {
    4398                   sprintf(generalPrint,"Testing OsiObject options %d\n",testOsiOptions);
     4976                  sprintf(generalPrint,"Testing OsiObject options %d",testOsiOptions);
    43994977                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
    44004978                    << generalPrint
     
    47115289                    babModel->makeGlobalCut(rowCutPointer);
    47125290                  }
    4713                 }
     5291                } 
    47145292                // If defaults then increase trust for small models
    47155293                if (!strongChanged) {
     
    47275305                babModel->setThreadMode(numberThreads/100);
    47285306#endif
     5307                int returnCode=callBack(babModel,3);
     5308                if (returnCode) {
     5309                  // exit if user wants
     5310                  delete babModel;
     5311                  return returnCode;
     5312                }
    47295313                babModel->branchAndBound(statistics);
     5314                returnCode=callBack(babModel,4);
     5315                if (returnCode) {
     5316                  // exit if user wants
     5317                  model.moveInfo(*babModel);
     5318                  delete babModel;
     5319                  return returnCode;
     5320                }
    47305321#ifdef CLP_MALLOC_STATISTICS
    47315322                malloc_stats();
     
    47665357                   could be passed to CbcClpUnitTest. */
    47675358                if (call_CbcClpUnitTest_on_777 == 777) {
    4768                   void CbcClpUnitTest (const CbcModel & saveModel,
    4769                                        std::string& dirMiplib);
    4770                   CbcClpUnitTest(model, dirMiplib);
     5359                  CbcClpUnitTest(model, dirMiplib, false);
    47715360                }
    47725361                return 777;
     
    48505439                }
    48515440                checkSOS(babModel, babModel->solver());
     5441              } else if (model.bestSolution()&&type==BAB&&model.getMinimizationObjValue()&&preProcess) {
     5442                sprintf(generalPrint,"Restoring heuristic best solution of %g",model.getMinimizationObjValue());
     5443                generalMessageHandler->message(CLP_GENERAL,generalMessages)
     5444                  << generalPrint
     5445                  <<CoinMessageEol;
     5446                int n = saveSolver->getNumCols();
     5447                bestSolution = new double [n];
     5448                // Put solution now back in saveSolver
     5449                babModel->assignSolver(saveSolver);
     5450                saveSolver->setColSolution(model.bestSolution());
     5451                memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
     5452                // and put back in very original solver
     5453                {
     5454                  ClpSimplex * original = originalSolver->getModelPtr();
     5455                  double * lower = original->columnLower();
     5456                  double * upper = original->columnUpper();
     5457                  double * solution = original->primalColumnSolution();
     5458                  int n = original->numberColumns();
     5459                  //assert (!n||n==babModel->solver()->getNumCols());
     5460                  for (int i=0;i<n;i++) {
     5461                    solution[i]=bestSolution[i];
     5462                    if (originalSolver->isInteger(i)) {
     5463                      lower[i]=solution[i];
     5464                      upper[i]=solution[i];
     5465                    }
     5466                  }
     5467                }
    48525468              }
    48535469              if (type==STRENGTHEN&&strengthenedModel)
     
    48875503                int iStat2 = babModel->secondaryStatus();
    48885504                if (!noPrinting) {
    4889                   sprintf(generalPrint,"Result - %s%s objective %.16g after %d nodes and %d iterations - took %.2f seconds",
     5505                  sprintf(generalPrint,"Result - %s%s objective %.16g after %d nodes and %d iterations - took %.2f seconds (total time %.2f)",
    48905506                          statusName[iStat].c_str(),minor[iStat2].c_str(),
    48915507                          babModel->getObjValue(),babModel->getNodeCount(),
    4892                           babModel->getIterationCount(),time2-time1);
     5508                          babModel->getIterationCount(),time2-time1,time2-time0);
    48935509                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
    48945510                    << generalPrint
    48955511                    <<CoinMessageEol;
     5512                }
     5513                int returnCode=callBack(babModel,5);
     5514                if (returnCode) {
     5515                  // exit if user wants
     5516                  model.moveInfo(*babModel);
     5517                  delete babModel;
     5518                  return returnCode;
    48965519                }
    48975520#ifdef COIN_HAS_ASL
     
    52845907                if (dualize) {
    52855908                  model2 = ((ClpSimplexOther *) model2)->dualOfModel();
    5286                   sprintf(generalPrint,"Dual of model has %d rows and %d columns\n",
     5909                  sprintf(generalPrint,"Dual of model has %d rows and %d columns",
    52875910                         model2->numberRows(),model2->numberColumns());
    52885911                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
     
    60586681          case UNITTEST:
    60596682            {
    6060               printf("unit test is now only from clp - does same thing\n");
    6061               //return(22);
     6683              CbcClpUnitTest(model, dirSample, true);
    60626684            }
    60636685            break;
     
    65677189  // By now all memory should be freed
    65687190#ifdef DMALLOC
    6569   dmalloc_log_unfreed();
    6570   dmalloc_shutdown();
     7191  //dmalloc_log_unfreed();
     7192  //dmalloc_shutdown();
    65717193#endif
    6572   if (babModel)
     7194  if (babModel) {
    65737195    model.moveInfo(*babModel);
    6574   delete babModel;
     7196    //babModel->setModelOwnsSolver(false);
     7197  }
     7198#ifdef CBC_SIG_TRAP
     7199  // On Sun sometimes seems to be error - try and get round it
     7200  CoinSighandler_t saveSignal=SIG_DFL;
     7201  // register signal handler
     7202  saveSignal=signal(SIGSEGV,signal_handler_error);
     7203  // to force failure!babModel->setNumberObjects(20000);
     7204  if (!sigsetjmp(cbc_seg_buffer,1)) {
     7205#endif
     7206    delete babModel;
     7207#ifdef CBC_SIG_TRAP
     7208  } else {
     7209    std::cerr<<"delete babModel failed"<<std::endl;
     7210  }
     7211  #endif
    65757212  model.solver()->setWarmStart(NULL);
     7213  sprintf(generalPrint,"Total time %.2f",CoinCpuTime()-time0);
     7214  generalMessageHandler->message(CLP_GENERAL,generalMessages)
     7215    << generalPrint
     7216    <<CoinMessageEol;
    65767217  return 0;
    65777218}   
  • branches/BSP/trunk/Cbc/src/unitTestClp.cpp

    r749 r764  
    22// Corporation and others.  All Rights Reserved.
    33
    4 #include "CoinPragma.hpp"
    5 #include <cassert>
    64#include <cstdio>
    7 #include <cmath>
    8 #include <cfloat>
    95#include <string>
    106#include <iostream>
    117
    12 #include "CoinMpsIO.hpp"
    13 #include "CoinPackedMatrix.hpp"
    14 #include "CoinPackedVector.hpp"
    15 #include "CoinHelperFunctions.hpp"
    168#include "CoinTime.hpp"
    179#include "CbcModel.hpp"
    1810#include "CbcCutGenerator.hpp"
    1911#include "OsiClpSolverInterface.hpp"
    20 #include "ClpFactorization.hpp"
    21 #include "ClpSimplex.hpp"
    22 #include "ClpSimplexOther.hpp"
    23 #include "ClpInterior.hpp"
    24 #include "ClpLinearObjective.hpp"
    25 #include "ClpDualRowSteepest.hpp"
    26 #include "ClpDualRowDantzig.hpp"
    27 #include "ClpPrimalColumnSteepest.hpp"
    28 #include "ClpPrimalColumnDantzig.hpp"
    29 #include "ClpParameters.hpp"
    30 #include "ClpNetworkMatrix.hpp"
    31 #include "ClpPlusMinusOneMatrix.hpp"
    32 #include "MyMessageHandler.hpp"
    33 #include "MyEventHandler.hpp"
    34 
    35 #include "ClpPresolve.hpp"
    36 #include "Idiot.hpp"
    37 
    3812
    3913//#############################################################################
     
    4317#endif
    4418
    45 // Function Prototypes. Function definitions is in this file.
    46 void testingMessage( const char * const msg );
    47 #if UFL_BARRIER
    48 static int barrierAvailable=1;
    49 static std::string nameBarrier="barrier-UFL";
    50 #elif WSSMP_BARRIER
    51 static int barrierAvailable=2;
    52 static std::string nameBarrier="barrier-WSSMP";
    53 #else
    54 static int barrierAvailable=0;
    55 static std::string nameBarrier="barrier-slow";
    56 #endif
    57 #define NUMBER_ALGORITHMS 12
    58 // If you just want a subset then set some to 1
    59 static int switchOff[NUMBER_ALGORITHMS]={0,0,0,0,0,0,0,0,0,0,0,0};
    60 // shortName - 0 no , 1 yes
    61 ClpSolve setupForSolve(int algorithm, std::string & nameAlgorithm,
    62                        int shortName)
    63 {
    64   ClpSolve solveOptions;
    65   /* algorithms are
    66      0 barrier
    67      1 dual with volumne crash
    68      2,3 dual with and without crash
    69      4,5 primal with and without
    70      6,7 automatic with and without
    71      8,9 primal with idiot 1 and 5
    72      10,11 primal with 70, dual with volume
    73   */
    74   switch (algorithm) {
    75   case 0:
    76     if (shortName)
    77       nameAlgorithm="ba";
    78     else
    79       nameAlgorithm="nameBarrier";
    80     solveOptions.setSolveType(ClpSolve::useBarrier);
    81     if (barrierAvailable==1)
    82       solveOptions.setSpecialOption(4,4);
    83     else if (barrierAvailable==2)
    84       solveOptions.setSpecialOption(4,2);
    85     break;
    86   case 1:
    87 #ifdef COIN_HAS_VOL
    88     if (shortName)
    89       nameAlgorithm="du-vol-50";
    90     else
    91       nameAlgorithm="dual-volume-50";
    92     solveOptions.setSolveType(ClpSolve::useDual);
    93     solveOptions.setSpecialOption(0,2,50); // volume
    94 #else 
    95       solveOptions.setSolveType(ClpSolve::notImplemented);
    96 #endif
    97     break;
    98   case 2:
    99     if (shortName)
    100       nameAlgorithm="du-cr";
    101     else
    102       nameAlgorithm="dual-crash";
    103     solveOptions.setSolveType(ClpSolve::useDual);
    104     solveOptions.setSpecialOption(0,1);
    105     break;
    106   case 3:
    107     if (shortName)
    108       nameAlgorithm="du";
    109     else
    110       nameAlgorithm="dual";
    111     solveOptions.setSolveType(ClpSolve::useDual);
    112     break;
    113   case 4:
    114     if (shortName)
    115       nameAlgorithm="pr-cr";
    116     else
    117       nameAlgorithm="primal-crash";
    118     solveOptions.setSolveType(ClpSolve::usePrimal);
    119     solveOptions.setSpecialOption(1,1);
    120     break;
    121   case 5:
    122     if (shortName)
    123       nameAlgorithm="pr";
    124     else
    125       nameAlgorithm="primal";
    126     solveOptions.setSolveType(ClpSolve::usePrimal);
    127     break;
    128   case 6:
    129     if (shortName)
    130       nameAlgorithm="au-cr";
    131     else
    132       nameAlgorithm="either-crash";
    133     solveOptions.setSolveType(ClpSolve::automatic);
    134     solveOptions.setSpecialOption(1,1);
    135     break;
    136   case 7:
    137     if (shortName)
    138       nameAlgorithm="au";
    139     else
    140       nameAlgorithm="either";
    141     solveOptions.setSolveType(ClpSolve::automatic);
    142     break;
    143   case 8:
    144     if (shortName)
    145       nameAlgorithm="pr-id-1";
    146     else
    147       nameAlgorithm="primal-idiot-1";
    148     solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    149     solveOptions.setSpecialOption(1,2,1); // idiot
    150     break;
    151   case 9:
    152     if (shortName)
    153       nameAlgorithm="pr-id-5";
    154     else
    155       nameAlgorithm="primal-idiot-5";
    156     solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    157     solveOptions.setSpecialOption(1,2,5); // idiot
    158     break;
    159   case 10:
    160     if (shortName)
    161       nameAlgorithm="pr-id-70";
    162     else
    163       nameAlgorithm="primal-idiot-70";
    164     solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    165     solveOptions.setSpecialOption(1,2,70); // idiot
    166     break;
    167   case 11:
    168 #ifdef COIN_HAS_VOL
    169     if (shortName)
    170       nameAlgorithm="du-vol";
    171     else
    172       nameAlgorithm="dual-volume";
    173     solveOptions.setSolveType(ClpSolve::useDual);
    174     solveOptions.setSpecialOption(0,2,3000); // volume
    175 #else 
    176       solveOptions.setSolveType(ClpSolve::notImplemented);
    177 #endif
    178     break;
    179   default:
    180     abort();
    181   }
    182   if (shortName) {
    183     // can switch off
    184     if (switchOff[algorithm])
    185       solveOptions.setSolveType(ClpSolve::notImplemented);
    186   }
    187   return solveOptions;
    188 }
    189 //----------------------------------------------------------------
    190 // unitTest [-mpsDir=V1] [-netlibDir=V2] [-test]
    191 //
    192 // where:
    193 //   -mpsDir: directory containing mps test files
    194 //       Default value V1="../../Data/Sample"   
    195 //   -netlibDir: directory containing netlib files
    196 //       Default value V2="../../Data/Netlib"
    197 //   -test
    198 //       If specified, then netlib test set run
    199 //
    200 // All parameters are optional.
    201 //----------------------------------------------------------------
    202 int mainTest (int argc, const char *argv[],int algorithm,
    203               ClpSimplex empty, bool doPresolve, int switchOffValue,bool doVector)
    204 {
    205   int i;
    206 
    207   if (switchOffValue>0) {
    208     // switch off some
    209     int iTest;
    210     for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
    211       int bottom = switchOffValue%10;
    212       assert (bottom==0||bottom==1);
    213       switchOffValue/= 10;
    214       switchOff[iTest]=0;
    215     }
    216   }
    217 
    218   // define valid parameter keywords
    219   std::set<std::string> definedKeyWords;
    220   definedKeyWords.insert("-mpsDir");
    221   definedKeyWords.insert("-netlibDir");
    222   definedKeyWords.insert("-netlib");
    223 
    224   // Create a map of parameter keys and associated data
    225   std::map<std::string,std::string> parms;
    226   for ( i=1; i<argc; i++ ) {
    227     std::string parm(argv[i]);
    228     std::string key,value;
    229     unsigned int  eqPos = parm.find('=');
    230 
    231     // Does parm contain and '='
    232     if ( eqPos==std::string::npos ) {
    233       //Parm does not contain '='
    234       key = parm;
    235     }
    236     else {
    237       key=parm.substr(0,eqPos);
    238       value=parm.substr(eqPos+1);
    239     }
    240 
    241     // Is specifed key valid?
    242     if ( definedKeyWords.find(key) == definedKeyWords.end() ) {
    243       // invalid key word.
    244       // Write help text
    245       std::cerr <<"Undefined parameter \"" <<key <<"\".\n";
    246       std::cerr <<"Correct usage: \n";
    247       std::cerr <<"  unitTest [-mpsDir=V1] [-netlibDir=V2] [-test[=V3]]\n";
    248       std::cerr <<"  where:\n";
    249       std::cerr <<"    -mpsDir: directory containing mps test files\n";
    250       std::cerr <<"        Default value V1=\"../../Data/Sample\"\n";
    251       std::cerr <<"    -netlibDir: directory containing netlib files\n";
    252       std::cerr <<"        Default value V2=\"../../Data/Netlib\"\n";
    253       std::cerr <<"    -test\n";
    254       std::cerr <<"        If specified, then netlib testset run.\n";
    255       std::cerr <<"        If V3 then taken as single file\n";
    256       return 1;
    257     }
    258     parms[key]=value;
    259   }
    260  
    261   const char dirsep =  CoinFindDirSeparator();
    262   // Set directory containing mps data files.
    263   std::string mpsDir;
    264   if (parms.find("-mpsDir") != parms.end())
    265     mpsDir=parms["-mpsDir"] + dirsep;
    266   else
    267     mpsDir = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\";
    268  
    269   // Set directory containing netlib data files.
    270   std::string netlibDir;
    271   if (parms.find("-netlibDir") != parms.end())
    272     netlibDir=parms["-netlibDir"] + dirsep;
    273   else
    274     netlibDir = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
    275   if (!empty.numberRows()) {
    276     testingMessage( "Testing ClpSimplex\n" );
    277     ClpSimplexUnitTest(mpsDir,netlibDir);
    278   }
    279   if (parms.find("-netlib") != parms.end()||empty.numberRows())
    280   {
    281     unsigned int m;
    282    
    283     // Define test problems:
    284     //   mps names,
    285     //   maximization or minimization,
    286     //   Number of rows and columns in problem, and
    287     //   objective function value
    288     std::vector<std::string> mpsName;
    289     std::vector<bool> min;
    290     std::vector<int> nRows;
    291     std::vector<int> nCols;
    292     std::vector<double> objValue;
    293     std::vector<double> objValueTol;
    294     // 100 added means no presolve
    295     std::vector<int> bestStrategy;
    296     if(empty.numberRows()) {
    297       std::string alg;
    298       for (int iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
    299         ClpSolve solveOptions=setupForSolve(iTest,alg,0);
    300         printf("%d %s ",iTest,alg.c_str());
    301         if (switchOff[iTest])
    302           printf("skipped by user\n");
    303         else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
    304           printf("skipped as not available\n");
    305         else
    306           printf("will be tested\n");
    307       }
    308     }
    309     if (!empty.numberRows()) {
    310       mpsName.push_back("25fv47");
    311       min.push_back(true);
    312       nRows.push_back(822);
    313       nCols.push_back(1571);
    314       objValueTol.push_back(1.E-10);
    315       objValue.push_back(5.5018458883E+03);
    316       bestStrategy.push_back(0);
    317       mpsName.push_back("80bau3b");min.push_back(true);nRows.push_back(2263);nCols.push_back(9799);objValueTol.push_back(1.e-10);objValue.push_back(9.8722419241E+05);bestStrategy.push_back(3);
    318       mpsName.push_back("blend");min.push_back(true);nRows.push_back(75);nCols.push_back(83);objValueTol.push_back(1.e-10);objValue.push_back(-3.0812149846e+01);bestStrategy.push_back(3);
    319       mpsName.push_back("pilotnov");min.push_back(true);nRows.push_back(976);nCols.push_back(2172);objValueTol.push_back(1.e-10);objValue.push_back(-4.4972761882e+03);bestStrategy.push_back(3);
    320       mpsName.push_back("maros-r7");min.push_back(true);nRows.push_back(3137);nCols.push_back(9408);objValueTol.push_back(1.e-10);objValue.push_back(1.4971851665e+06);bestStrategy.push_back(2);
    321       mpsName.push_back("pilot");min.push_back(true);nRows.push_back(1442);nCols.push_back(3652);objValueTol.push_back(1.e-5);objValue.push_back(/*-5.5740430007e+02*/-557.48972927292);bestStrategy.push_back(3);
    322       mpsName.push_back("pilot4");min.push_back(true);nRows.push_back(411);nCols.push_back(1000);objValueTol.push_back(5.e-5);objValue.push_back(-2.5811392641e+03);bestStrategy.push_back(3);
    323       mpsName.push_back("pilot87");min.push_back(true);nRows.push_back(2031);nCols.push_back(4883);objValueTol.push_back(1.e-4);objValue.push_back(3.0171072827e+02);bestStrategy.push_back(0);
    324       mpsName.push_back("adlittle");min.push_back(true);nRows.push_back(57);nCols.push_back(97);objValueTol.push_back(1.e-10);objValue.push_back(2.2549496316e+05);bestStrategy.push_back(3);
    325       mpsName.push_back("afiro");min.push_back(true);nRows.push_back(28);nCols.push_back(32);objValueTol.push_back(1.e-10);objValue.push_back(-4.6475314286e+02);bestStrategy.push_back(3);
    326       mpsName.push_back("agg");min.push_back(true);nRows.push_back(489);nCols.push_back(163);objValueTol.push_back(1.e-10);objValue.push_back(-3.5991767287e+07);bestStrategy.push_back(3);
    327       mpsName.push_back("agg2");min.push_back(true);nRows.push_back(517);nCols.push_back(302);objValueTol.push_back(1.e-10);objValue.push_back(-2.0239252356e+07);bestStrategy.push_back(3);
    328       mpsName.push_back("agg3");min.push_back(true);nRows.push_back(517);nCols.push_back(302);objValueTol.push_back(1.e-10);objValue.push_back(1.0312115935e+07);bestStrategy.push_back(4);
    329       mpsName.push_back("bandm");min.push_back(true);nRows.push_back(306);nCols.push_back(472);objValueTol.push_back(1.e-10);objValue.push_back(-1.5862801845e+02);bestStrategy.push_back(2);
    330       mpsName.push_back("beaconfd");min.push_back(true);nRows.push_back(174);nCols.push_back(262);objValueTol.push_back(1.e-10);objValue.push_back(3.3592485807e+04);bestStrategy.push_back(0);
    331       mpsName.push_back("bnl1");min.push_back(true);nRows.push_back(644);nCols.push_back(1175);objValueTol.push_back(1.e-10);objValue.push_back(1.9776295615E+03);bestStrategy.push_back(3);
    332       mpsName.push_back("bnl2");min.push_back(true);nRows.push_back(2325);nCols.push_back(3489);objValueTol.push_back(1.e-10);objValue.push_back(1.8112365404e+03);bestStrategy.push_back(3);
    333       mpsName.push_back("boeing1");min.push_back(true);nRows.push_back(/*351*/352);nCols.push_back(384);objValueTol.push_back(1.e-10);objValue.push_back(-3.3521356751e+02);bestStrategy.push_back(3);
    334       mpsName.push_back("boeing2");min.push_back(true);nRows.push_back(167);nCols.push_back(143);objValueTol.push_back(1.e-10);objValue.push_back(-3.1501872802e+02);bestStrategy.push_back(3);
    335       mpsName.push_back("bore3d");min.push_back(true);nRows.push_back(234);nCols.push_back(315);objValueTol.push_back(1.e-10);objValue.push_back(1.3730803942e+03);bestStrategy.push_back(3);
    336       mpsName.push_back("brandy");min.push_back(true);nRows.push_back(221);nCols.push_back(249);objValueTol.push_back(1.e-10);objValue.push_back(1.5185098965e+03);bestStrategy.push_back(3);
    337       mpsName.push_back("capri");min.push_back(true);nRows.push_back(272);nCols.push_back(353);objValueTol.push_back(1.e-10);objValue.push_back(2.6900129138e+03);bestStrategy.push_back(3);
    338       mpsName.push_back("cycle");min.push_back(true);nRows.push_back(1904);nCols.push_back(2857);objValueTol.push_back(1.e-9);objValue.push_back(-5.2263930249e+00);bestStrategy.push_back(3);
    339       mpsName.push_back("czprob");min.push_back(true);nRows.push_back(930);nCols.push_back(3523);objValueTol.push_back(1.e-10);objValue.push_back(2.1851966989e+06);bestStrategy.push_back(3);
    340       mpsName.push_back("d2q06c");min.push_back(true);nRows.push_back(2172);nCols.push_back(5167);objValueTol.push_back(1.e-7);objValue.push_back(122784.21557456);bestStrategy.push_back(0);
    341       mpsName.push_back("d6cube");min.push_back(true);nRows.push_back(416);nCols.push_back(6184);objValueTol.push_back(1.e-7);objValue.push_back(3.1549166667e+02);bestStrategy.push_back(3);
    342       mpsName.push_back("degen2");min.push_back(true);nRows.push_back(445);nCols.push_back(534);objValueTol.push_back(1.e-10);objValue.push_back(-1.4351780000e+03);bestStrategy.push_back(3);
    343       mpsName.push_back("degen3");min.push_back(true);nRows.push_back(1504);nCols.push_back(1818);objValueTol.push_back(1.e-10);objValue.push_back(-9.8729400000e+02);bestStrategy.push_back(2);
    344       mpsName.push_back("dfl001");min.push_back(true);nRows.push_back(6072);nCols.push_back(12230);objValueTol.push_back(1.e-5);objValue.push_back(1.1266396047E+07);bestStrategy.push_back(5);
    345       mpsName.push_back("e226");min.push_back(true);nRows.push_back(224);nCols.push_back(282);objValueTol.push_back(1.e-10);objValue.push_back(-1.8751929066e+01+7.113);bestStrategy.push_back(3); // The correct answer includes -7.113 term. This is a constant in the objective function. See line 1683 of the mps file.
    346       mpsName.push_back("etamacro");min.push_back(true);nRows.push_back(401);nCols.push_back(688);objValueTol.push_back(1.e-6);objValue.push_back(-7.5571521774e+02 );bestStrategy.push_back(3);
    347       mpsName.push_back("fffff800");min.push_back(true);nRows.push_back(525);nCols.push_back(854);objValueTol.push_back(1.e-6);objValue.push_back(5.5567961165e+05);bestStrategy.push_back(3);
    348       mpsName.push_back("finnis");min.push_back(true);nRows.push_back(498);nCols.push_back(614);objValueTol.push_back(1.e-6);objValue.push_back(1.7279096547e+05);bestStrategy.push_back(3);
    349       mpsName.push_back("fit1d");min.push_back(true);nRows.push_back(25);nCols.push_back(1026);objValueTol.push_back(1.e-10);objValue.push_back(-9.1463780924e+03);bestStrategy.push_back(3+100);
    350       mpsName.push_back("fit1p");min.push_back(true);nRows.push_back(628);nCols.push_back(1677);objValueTol.push_back(1.e-10);objValue.push_back(9.1463780924e+03);bestStrategy.push_back(5+100);
    351       mpsName.push_back("fit2d");min.push_back(true);nRows.push_back(26);nCols.push_back(10500);objValueTol.push_back(1.e-10);objValue.push_back(-6.8464293294e+04);bestStrategy.push_back(3+100);
    352       mpsName.push_back("fit2p");min.push_back(true);nRows.push_back(3001);nCols.push_back(13525);objValueTol.push_back(1.e-9);objValue.push_back(6.8464293232e+04);bestStrategy.push_back(5+100);
    353       mpsName.push_back("forplan");min.push_back(true);nRows.push_back(162);nCols.push_back(421);objValueTol.push_back(1.e-6);objValue.push_back(-6.6421873953e+02);bestStrategy.push_back(3);
    354       mpsName.push_back("ganges");min.push_back(true);nRows.push_back(1310);nCols.push_back(1681);objValueTol.push_back(1.e-5);objValue.push_back(-1.0958636356e+05);bestStrategy.push_back(3);
    355       mpsName.push_back("gfrd-pnc");min.push_back(true);nRows.push_back(617);nCols.push_back(1092);objValueTol.push_back(1.e-10);objValue.push_back(6.9022359995e+06);bestStrategy.push_back(3);
    356       mpsName.push_back("greenbea");min.push_back(true);nRows.push_back(2393);nCols.push_back(5405);objValueTol.push_back(1.e-10);objValue.push_back(/*-7.2462405908e+07*/-72555248.129846);bestStrategy.push_back(3);
    357       mpsName.push_back("greenbeb");min.push_back(true);nRows.push_back(2393);nCols.push_back(5405);objValueTol.push_back(1.e-10);objValue.push_back(/*-4.3021476065e+06*/-4302260.2612066);bestStrategy.push_back(3);
    358       mpsName.push_back("grow15");min.push_back(true);nRows.push_back(301);nCols.push_back(645);objValueTol.push_back(1.e-10);objValue.push_back(-1.0687094129e+08);bestStrategy.push_back(4+100);
    359       mpsName.push_back("grow22");min.push_back(true);nRows.push_back(441);nCols.push_back(946);objValueTol.push_back(1.e-10);objValue.push_back(-1.6083433648e+08);bestStrategy.push_back(4+100);
    360       mpsName.push_back("grow7");min.push_back(true);nRows.push_back(141);nCols.push_back(301);objValueTol.push_back(1.e-10);objValue.push_back(-4.7787811815e+07);bestStrategy.push_back(4+100);
    361       mpsName.push_back("israel");min.push_back(true);nRows.push_back(175);nCols.push_back(142);objValueTol.push_back(1.e-10);objValue.push_back(-8.9664482186e+05);bestStrategy.push_back(2);
    362       mpsName.push_back("kb2");min.push_back(true);nRows.push_back(44);nCols.push_back(41);objValueTol.push_back(1.e-10);objValue.push_back(-1.7499001299e+03);bestStrategy.push_back(3);
    363       mpsName.push_back("lotfi");min.push_back(true);nRows.push_back(154);nCols.push_back(308);objValueTol.push_back(1.e-10);objValue.push_back(-2.5264706062e+01);bestStrategy.push_back(3);
    364       mpsName.push_back("maros");min.push_back(true);nRows.push_back(847);nCols.push_back(1443);objValueTol.push_back(1.e-10);objValue.push_back(-5.8063743701e+04);bestStrategy.push_back(3);
    365       mpsName.push_back("modszk1");min.push_back(true);nRows.push_back(688);nCols.push_back(1620);objValueTol.push_back(1.e-10);objValue.push_back(3.2061972906e+02);bestStrategy.push_back(3);
    366       mpsName.push_back("nesm");min.push_back(true);nRows.push_back(663);nCols.push_back(2923);objValueTol.push_back(1.e-5);objValue.push_back(1.4076073035e+07);bestStrategy.push_back(2);
    367       mpsName.push_back("perold");min.push_back(true);nRows.push_back(626);nCols.push_back(1376);objValueTol.push_back(1.e-6);objValue.push_back(-9.3807580773e+03);bestStrategy.push_back(3);
    368       //mpsName.push_back("qap12");min.push_back(true);nRows.push_back(3193);nCols.push_back(8856);objValueTol.push_back(1.e-6);objValue.push_back(5.2289435056e+02);bestStrategy.push_back(3);
    369       //mpsName.push_back("qap15");min.push_back(true);nRows.push_back(6331);nCols.push_back(22275);objValueTol.push_back(1.e-10);objValue.push_back(1.0409940410e+03);bestStrategy.push_back(3);
    370       mpsName.push_back("recipe");min.push_back(true);nRows.push_back(92);nCols.push_back(180);objValueTol.push_back(1.e-10);objValue.push_back(-2.6661600000e+02);bestStrategy.push_back(3);
    371       mpsName.push_back("sc105");min.push_back(true);nRows.push_back(106);nCols.push_back(103);objValueTol.push_back(1.e-10);objValue.push_back(-5.2202061212e+01);bestStrategy.push_back(3);
    372       mpsName.push_back("sc205");min.push_back(true);nRows.push_back(206);nCols.push_back(203);objValueTol.push_back(1.e-10);objValue.push_back(-5.2202061212e+01);bestStrategy.push_back(3);
    373       mpsName.push_back("sc50a");min.push_back(true);nRows.push_back(51);nCols.push_back(48);objValueTol.push_back(1.e-10);objValue.push_back(-6.4575077059e+01);bestStrategy.push_back(3);
    374       mpsName.push_back("sc50b");min.push_back(true);nRows.push_back(51);nCols.push_back(48);objValueTol.push_back(1.e-10);objValue.push_back(-7.0000000000e+01);bestStrategy.push_back(3);
    375       mpsName.push_back("scagr25");min.push_back(true);nRows.push_back(472);nCols.push_back(500);objValueTol.push_back(1.e-10);objValue.push_back(-1.4753433061e+07);bestStrategy.push_back(3);
    376       mpsName.push_back("scagr7");min.push_back(true);nRows.push_back(130);nCols.push_back(140);objValueTol.push_back(1.e-6);objValue.push_back(-2.3313892548e+06);bestStrategy.push_back(3);
    377       mpsName.push_back("scfxm1");min.push_back(true);nRows.push_back(331);nCols.push_back(457);objValueTol.push_back(1.e-10);objValue.push_back(1.8416759028e+04);bestStrategy.push_back(3);
    378       mpsName.push_back("scfxm2");min.push_back(true);nRows.push_back(661);nCols.push_back(914);objValueTol.push_back(1.e-10);objValue.push_back(3.6660261565e+04);bestStrategy.push_back(3);
    379       mpsName.push_back("scfxm3");min.push_back(true);nRows.push_back(991);nCols.push_back(1371);objValueTol.push_back(1.e-10);objValue.push_back(5.4901254550e+04);bestStrategy.push_back(3);
    380       mpsName.push_back("scorpion");min.push_back(true);nRows.push_back(389);nCols.push_back(358);objValueTol.push_back(1.e-10);objValue.push_back(1.8781248227e+03);bestStrategy.push_back(3);
    381       mpsName.push_back("scrs8");min.push_back(true);nRows.push_back(491);nCols.push_back(1169);objValueTol.push_back(1.e-5);objValue.push_back(9.0429998619e+02);bestStrategy.push_back(2);
    382       mpsName.push_back("scsd1");min.push_back(true);nRows.push_back(78);nCols.push_back(760);objValueTol.push_back(1.e-10);objValue.push_back(8.6666666743e+00);bestStrategy.push_back(3+100);
    383       mpsName.push_back("scsd6");min.push_back(true);nRows.push_back(148);nCols.push_back(1350);objValueTol.push_back(1.e-10);objValue.push_back(5.0500000078e+01);bestStrategy.push_back(3+100);
    384       mpsName.push_back("scsd8");min.push_back(true);nRows.push_back(398);nCols.push_back(2750);objValueTol.push_back(1.e-10);objValue.push_back(9.0499999993e+02);bestStrategy.push_back(1+100);
    385       mpsName.push_back("sctap1");min.push_back(true);nRows.push_back(301);nCols.push_back(480);objValueTol.push_back(1.e-10);objValue.push_back(1.4122500000e+03);bestStrategy.push_back(3);
    386       mpsName.push_back("sctap2");min.push_back(true);nRows.push_back(1091);nCols.push_back(1880);objValueTol.push_back(1.e-10);objValue.push_back(1.7248071429e+03);bestStrategy.push_back(3);
    387       mpsName.push_back("sctap3");min.push_back(true);nRows.push_back(1481);nCols.push_back(2480);objValueTol.push_back(1.e-10);objValue.push_back(1.4240000000e+03);bestStrategy.push_back(3);
    388       mpsName.push_back("seba");min.push_back(true);nRows.push_back(516);nCols.push_back(1028);objValueTol.push_back(1.e-10);objValue.push_back(1.5711600000e+04);bestStrategy.push_back(3);
    389       mpsName.push_back("share1b");min.push_back(true);nRows.push_back(118);nCols.push_back(225);objValueTol.push_back(1.e-10);objValue.push_back(-7.6589318579e+04);bestStrategy.push_back(3);
    390       mpsName.push_back("share2b");min.push_back(true);nRows.push_back(97);nCols.push_back(79);objValueTol.push_back(1.e-10);objValue.push_back(-4.1573224074e+02);bestStrategy.push_back(3);
    391       mpsName.push_back("shell");min.push_back(true);nRows.push_back(537);nCols.push_back(1775);objValueTol.push_back(1.e-10);objValue.push_back(1.2088253460e+09);bestStrategy.push_back(3);
    392       mpsName.push_back("ship04l");min.push_back(true);nRows.push_back(403);nCols.push_back(2118);objValueTol.push_back(1.e-10);objValue.push_back(1.7933245380e+06);bestStrategy.push_back(3);
    393       mpsName.push_back("ship04s");min.push_back(true);nRows.push_back(403);nCols.push_back(1458);objValueTol.push_back(1.e-10);objValue.push_back(1.7987147004e+06);bestStrategy.push_back(3);
    394       mpsName.push_back("ship08l");min.push_back(true);nRows.push_back(779);nCols.push_back(4283);objValueTol.push_back(1.e-10);objValue.push_back(1.9090552114e+06);bestStrategy.push_back(3);
    395       mpsName.push_back("ship08s");min.push_back(true);nRows.push_back(779);nCols.push_back(2387);objValueTol.push_back(1.e-10);objValue.push_back(1.9200982105e+06);bestStrategy.push_back(2);
    396       mpsName.push_back("ship12l");min.push_back(true);nRows.push_back(1152);nCols.push_back(5427);objValueTol.push_back(1.e-10);objValue.push_back(1.4701879193e+06);bestStrategy.push_back(3);
    397       mpsName.push_back("ship12s");min.push_back(true);nRows.push_back(1152);nCols.push_back(2763);objValueTol.push_back(1.e-10);objValue.push_back(1.4892361344e+06);bestStrategy.push_back(2);
    398       mpsName.push_back("sierra");min.push_back(true);nRows.push_back(1228);nCols.push_back(2036);objValueTol.push_back(1.e-10);objValue.push_back(1.5394362184e+07);bestStrategy.push_back(3);
    399       mpsName.push_back("stair");min.push_back(true);nRows.push_back(357);nCols.push_back(467);objValueTol.push_back(1.e-10);objValue.push_back(-2.5126695119e+02);bestStrategy.push_back(3);
    400       mpsName.push_back("standata");min.push_back(true);nRows.push_back(360);nCols.push_back(1075);objValueTol.push_back(1.e-10);objValue.push_back(1.2576995000e+03);bestStrategy.push_back(3);
    401       //mpsName.push_back("standgub");min.push_back(true);nRows.push_back(362);nCols.push_back(1184);objValueTol.push_back(1.e-10);objValue.push_back(1257.6995); bestStrategy.push_back(3);
    402       mpsName.push_back("standmps");min.push_back(true);nRows.push_back(468);nCols.push_back(1075);objValueTol.push_back(1.e-10);objValue.push_back(1.4060175000E+03); bestStrategy.push_back(3);
    403       mpsName.push_back("stocfor1");min.push_back(true);nRows.push_back(118);nCols.push_back(111);objValueTol.push_back(1.e-10);objValue.push_back(-4.1131976219E+04);bestStrategy.push_back(3);
    404       mpsName.push_back("stocfor2");min.push_back(true);nRows.push_back(2158);nCols.push_back(2031);objValueTol.push_back(1.e-10);objValue.push_back(-3.9024408538e+04);bestStrategy.push_back(3);
    405       //mpsName.push_back("stocfor3");min.push_back(true);nRows.push_back(16676);nCols.push_back(15695);objValueTol.push_back(1.e-10);objValue.push_back(-3.9976661576e+04);bestStrategy.push_back(3);
    406       //mpsName.push_back("truss");min.push_back(true);nRows.push_back(1001);nCols.push_back(8806);objValueTol.push_back(1.e-10);objValue.push_back(4.5881584719e+05);bestStrategy.push_back(3);
    407       mpsName.push_back("tuff");min.push_back(true);nRows.push_back(334);nCols.push_back(587);objValueTol.push_back(1.e-10);objValue.push_back(2.9214776509e-01);bestStrategy.push_back(3);
    408       mpsName.push_back("vtpbase");min.push_back(true);nRows.push_back(199);nCols.push_back(203);objValueTol.push_back(1.e-10);objValue.push_back(1.2983146246e+05);bestStrategy.push_back(3);
    409       mpsName.push_back("wood1p");min.push_back(true);nRows.push_back(245);nCols.push_back(2594);objValueTol.push_back(5.e-5);objValue.push_back(1.4429024116e+00);bestStrategy.push_back(3);
    410       mpsName.push_back("woodw");min.push_back(true);nRows.push_back(1099);nCols.push_back(8405);objValueTol.push_back(1.e-10);objValue.push_back(1.3044763331E+00);bestStrategy.push_back(3);
    411     } else {
    412       // Just testing one
    413       mpsName.push_back(empty.problemName());min.push_back(true);nRows.push_back(-1);
    414       nCols.push_back(-1);objValueTol.push_back(1.e-10);
    415       objValue.push_back(0.0);bestStrategy.push_back(0);
    416       int iTest;
    417       std::string alg;
    418       for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
    419         ClpSolve solveOptions=setupForSolve(iTest,alg,0);
    420         printf("%d %s ",iTest,alg.c_str());
    421         if (switchOff[iTest])
    422           printf("skipped by user\n");
    423         else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
    424           printf("skipped as not available\n");
    425         else
    426           printf("will be tested\n");
    427       }
    428     }
    429 
    430     double timeTaken =0.0;
    431     if( !barrierAvailable)
    432       switchOff[0]=1;
    433   // Loop once for each Mps File
    434     for (m=0; m<mpsName.size(); m++ ) {
    435       std::cerr <<"  processing mps file: " <<mpsName[m]
    436                 <<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl;
    437 
    438       ClpSimplex solutionBase=empty;
    439       std::string fn = netlibDir+mpsName[m];
    440       if (!empty.numberRows()||algorithm<6) {
    441         // Read data mps file,
    442         CoinMpsIO mps;
    443         mps.readMps(fn.c_str(),"mps");
    444         solutionBase.loadProblem(*mps.getMatrixByCol(),mps.getColLower(),
    445                                  mps.getColUpper(),
    446                                  mps.getObjCoefficients(),
    447                                  mps.getRowLower(),mps.getRowUpper());
    448        
    449         solutionBase.setDblParam(ClpObjOffset,mps.objectiveOffset());
    450       }
    451      
    452       // Runs through strategies
    453       if (algorithm==6||algorithm==7) {
    454         // algorithms tested are at top of file
    455         double testTime[NUMBER_ALGORITHMS];
    456         std::string alg[NUMBER_ALGORITHMS];
    457         int iTest;
    458         for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
    459           ClpSolve solveOptions=setupForSolve(iTest,alg[iTest],1);
    460           if (solveOptions.getSolveType()!=ClpSolve::notImplemented) {
    461             double time1 = CoinCpuTime();
    462             ClpSimplex solution=solutionBase;
    463             if (solution.maximumSeconds()<0.0)
    464               solution.setMaximumSeconds(120.0);
    465             if (doVector) {
    466               ClpMatrixBase * matrix = solution.clpMatrix();
    467               if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
    468                 ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    469                 clpMatrix->makeSpecialColumnCopy();
    470               }
    471             }
    472             solution.initialSolve(solveOptions);
    473             double time2 = CoinCpuTime()-time1;
    474             testTime[iTest]=time2;
    475             printf("Took %g seconds - status %d\n",time2,solution.problemStatus());
    476             if (solution.problemStatus())
    477               testTime[iTest]=1.0e20;
    478           } else {
    479             testTime[iTest]=1.0e30;
    480           }
    481         }
    482         int iBest=-1;
    483         double dBest=1.0e10;
    484         printf("%s",fn.c_str());
    485         for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
    486           if (testTime[iTest]<1.0e30) {
    487             printf(" %s %g",
    488                    alg[iTest].c_str(),testTime[iTest]);
    489             if (testTime[iTest]<dBest) {
    490               dBest=testTime[iTest];
    491               iBest=iTest;
    492             }
    493           }
    494         }
    495         printf("\n");
    496         if (iBest>=0)
    497           printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
    498                  fn.c_str(),alg[iBest].c_str(),iBest,testTime[iBest]);
    499         else
    500           printf("No strategy finished in time limit\n");
    501         continue;
    502       }
    503       double time1 = CoinCpuTime();
    504       ClpSimplex solution=solutionBase;
    505 #if 0
    506       solution.setOptimizationDirection(-1);
    507       {
    508         int j;
    509         double * obj = solution.objective();
    510         int n=solution.numberColumns();
    511         for (j=0;j<n;j++)
    512           obj[j] *= -1.0;
    513       }
    514 #endif
    515       ClpSolve::SolveType method;
    516       ClpSolve::PresolveType presolveType;
    517       ClpSolve solveOptions;
    518       if (doPresolve)
    519         presolveType=ClpSolve::presolveOn;
    520       else
    521         presolveType=ClpSolve::presolveOff;
    522       solveOptions.setPresolveType(presolveType,5);
    523       std::string nameAlgorithm;
    524       if (algorithm!=5) {
    525         if (algorithm==0) {
    526           method=ClpSolve::useDual;
    527           nameAlgorithm="dual";
    528         } else if (algorithm==1) {
    529           method=ClpSolve::usePrimalorSprint;
    530           nameAlgorithm="primal";
    531         } else if (algorithm==3) {
    532           method=ClpSolve::automatic;
    533           nameAlgorithm="either";
    534         } else {
    535           nameAlgorithm="barrier-slow";
    536 #ifdef UFL_BARRIER
    537           solveOptions.setSpecialOption(4,4);
    538           nameAlgorithm="barrier-UFL";
    539 #endif
    540 #ifdef WSSMP_BARRIER
    541           solveOptions.setSpecialOption(4,2);
    542           nameAlgorithm="barrier-WSSMP";
    543 #endif
    544           method = ClpSolve::useBarrier;
    545         }
    546         solveOptions.setSolveType(method);
    547       } else {
    548         int iAlg = bestStrategy[m];
    549         int presolveOff=iAlg/100;
    550         iAlg=iAlg % 100;
    551         if( !barrierAvailable&&iAlg==0) {
    552           if (nRows[m]!=2172)
    553             iAlg = 5; // try primal
    554           else
    555             iAlg=3; // d2q06c
    556         }
    557         solveOptions=setupForSolve(iAlg,nameAlgorithm,0);
    558         if (presolveOff)
    559           solveOptions.setPresolveType(ClpSolve::presolveOff);
    560       }
    561       if (doVector) {
    562         ClpMatrixBase * matrix = solution.clpMatrix();
    563         if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
    564           ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    565           clpMatrix->makeSpecialColumnCopy();
    566         }
    567       }
    568       solution.initialSolve(solveOptions);
    569       double time2 = CoinCpuTime()-time1;
    570       timeTaken += time2;
    571       printf("%s took %g seconds using algorithm %s\n",fn.c_str(),time2,nameAlgorithm.c_str());
    572       // Test objective solution value
    573       {
    574         double soln = solution.objectiveValue();
    575         CoinRelFltEq eq(objValueTol[m]);
    576         std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
    577           soln-objValue[m]<<std::endl;
    578         if(!eq(soln,objValue[m]))
    579           printf("** difference fails\n");
    580       }
    581     }
    582     printf("Total time %g seconds\n",timeTaken);
    583   }
    584   else {
    585     testingMessage( "***Skipped Testing on netlib    ***\n" );
    586     testingMessage( "***use -netlib to test class***\n" );
    587   }
    588  
    589   testingMessage( "All tests completed successfully\n" );
    590   return 0;
    591 }
    592 
    593  
     19//#############################################################################
     20
    59421// Display message on stdout and stderr
    59522void testingMessage( const char * const msg )
    59623{
    597   std::cerr <<msg;
     24  std::cout <<msg;
    59825  //cout <<endl <<"*****************************************"
    59926  //     <<endl <<msg <<endl;
    60027}
    60128
    602 //--------------------------------------------------------------------------
    603 // test factorization methods and simplex method and simple barrier
    604 void
    605 ClpSimplexUnitTest(const std::string & mpsDir,
    606                    const std::string & netlibDir)
    607 {
    608  
    609   CoinRelFltEq eq(0.000001);
    610 
    611   {
    612     ClpSimplex solution;
    613  
    614     // matrix data
    615     //deliberate hiccup of 2 between 0 and 1
    616     CoinBigIndex start[5]={0,4,7,8,9};
    617     int length[5]={2,3,1,1,1};
    618     int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
    619     double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
    620     CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
    621    
    622     // rim data
    623     double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
    624     double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
    625     double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
    626     double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
    627     double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
    628    
    629     // basis 1
    630     int rowBasis1[3]={-1,-1,-1};
    631     int colBasis1[5]={1,1,-1,-1,1};
    632     solution.loadProblem(matrix,colLower,colUpper,objective,
    633                          rowLower,rowUpper);
    634     int i;
    635     solution.createStatus();
    636     for (i=0;i<3;i++) {
    637       if (rowBasis1[i]<0) {
    638         solution.setRowStatus(i,ClpSimplex::atLowerBound);
    639       } else {
    640         solution.setRowStatus(i,ClpSimplex::basic);
    641       }
    642     }
    643     for (i=0;i<5;i++) {
    644       if (colBasis1[i]<0) {
    645         solution.setColumnStatus(i,ClpSimplex::atLowerBound);
    646       } else {
    647         solution.setColumnStatus(i,ClpSimplex::basic);
    648       }
    649     }
    650     solution.setLogLevel(3+4+8+16+32);
    651     solution.primal();
    652     for (i=0;i<3;i++) {
    653       if (rowBasis1[i]<0) {
    654         solution.setRowStatus(i,ClpSimplex::atLowerBound);
    655       } else {
    656         solution.setRowStatus(i,ClpSimplex::basic);
    657       }
    658     }
    659     for (i=0;i<5;i++) {
    660       if (colBasis1[i]<0) {
    661         solution.setColumnStatus(i,ClpSimplex::atLowerBound);
    662       } else {
    663         solution.setColumnStatus(i,ClpSimplex::basic);
    664       }
    665     }
    666     // intricate stuff does not work with scaling
    667     solution.scaling(0);
    668     int returnCode = solution.factorize ( );
    669     assert(!returnCode);
    670     const double * colsol = solution.primalColumnSolution();
    671     const double * rowsol = solution.primalRowSolution();
    672     solution.getSolution(rowsol,colsol);
    673     double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
    674     for (i=0;i<5;i++) {
    675       assert(eq(colsol[i],colsol1[i]));
    676     }
    677     // now feed in again without actually doing factorization
    678     ClpFactorization factorization2 = *solution.factorization();
    679     ClpSimplex solution2 = solution;
    680     solution2.setFactorization(factorization2);
    681     solution2.createStatus();
    682     for (i=0;i<3;i++) {
    683       if (rowBasis1[i]<0) {
    684         solution2.setRowStatus(i,ClpSimplex::atLowerBound);
    685       } else {
    686         solution2.setRowStatus(i,ClpSimplex::basic);
    687       }
    688     }
    689     for (i=0;i<5;i++) {
    690       if (colBasis1[i]<0) {
    691         solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
    692       } else {
    693         solution2.setColumnStatus(i,ClpSimplex::basic);
    694       }
    695     }
    696     // intricate stuff does not work with scaling
    697     solution2.scaling(0);
    698     solution2.getSolution(rowsol,colsol);
    699     colsol = solution2.primalColumnSolution();
    700     rowsol = solution2.primalRowSolution();
    701     for (i=0;i<5;i++) {
    702       assert(eq(colsol[i],colsol1[i]));
    703     }
    704     solution2.setDualBound(0.1);
    705     solution2.dual();
    706     objective[2]=-1.0;
    707     objective[3]=-0.5;
    708     objective[4]=10.0;
    709     solution.dual();
    710     for (i=0;i<3;i++) {
    711       rowLower[i]=-1.0e20;
    712       colUpper[i+2]=0.0;
    713     }
    714     solution.setLogLevel(3);
    715     solution.dual();
    716     double rowObjective[]={1.0,0.5,-10.0};
    717     solution.loadProblem(matrix,colLower,colUpper,objective,
    718                          rowLower,rowUpper,rowObjective);
    719     solution.dual();
    720     solution.loadProblem(matrix,colLower,colUpper,objective,
    721                          rowLower,rowUpper,rowObjective);
    722     solution.primal();
    723   }
    724 #ifndef COIN_NO_CLP_MESSAGE
    725   {   
    726     CoinMpsIO m;
    727     std::string fn = mpsDir+"exmip1";
    728     m.readMps(fn.c_str(),"mps");
    729     ClpSimplex solution;
    730     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    731                          m.getObjCoefficients(),
    732                          m.getRowLower(),m.getRowUpper());
    733     solution.dual();
    734     // Test event handling
    735     MyEventHandler handler;
    736     solution.passInEventHandler(&handler);
    737     int numberRows=solution.numberRows();
    738     // make sure values pass has something to do
    739     for (int i=0;i<numberRows;i++)
    740       solution.setRowStatus(i,ClpSimplex::basic);
    741     solution.primal(1);
    742     assert (solution.secondaryStatus()==102); // Came out at end of pass
    743   }
    744   // Test Message handler
    745   {   
    746     CoinMpsIO m;
    747     std::string fn = mpsDir+"exmip1";
    748     //fn = "Test/subGams4";
    749     m.readMps(fn.c_str(),"mps");
    750     ClpSimplex model;
    751     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    752                          m.getObjCoefficients(),
    753                          m.getRowLower(),m.getRowUpper());
    754     // Message handler
    755     MyMessageHandler messageHandler(&model);
    756     std::cout<<"Testing derived message handler"<<std::endl;
    757     model.passInMessageHandler(&messageHandler);
    758     model.messagesPointer()->setDetailMessage(1,102);
    759     model.setFactorizationFrequency(10);
    760     model.primal();
    761 
    762     // Write saved solutions
    763     int nc = model.getNumCols();
    764     int s;
    765     std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
    766     int numSavedSolutions = fep.size();
    767     for ( s=0; s<numSavedSolutions; ++s ) {
    768       const StdVectorDouble & solnVec = fep[s];
    769       for ( int c=0; c<nc; ++c ) {
    770         if (fabs(solnVec[c])>1.0e-8)
    771           std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
    772       }
    773     }
    774     // Solve again without scaling
    775     // and maximize then minimize
    776     messageHandler.clearFeasibleExtremePoints();
    777     model.scaling(0);
    778     model.setOptimizationDirection(-1);
    779     model.primal();
    780     model.setOptimizationDirection(1);
    781     model.primal();
    782     fep = messageHandler.getFeasibleExtremePoints();
    783     numSavedSolutions = fep.size();
    784     for ( s=0; s<numSavedSolutions; ++s ) {
    785       const StdVectorDouble & solnVec = fep[s];
    786       for ( int c=0; c<nc; ++c ) {
    787         if (fabs(solnVec[c])>1.0e-8)
    788           std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
    789       }
    790     }
    791   }
    792 #endif
    793   // Test dual ranging
    794   {   
    795     CoinMpsIO m;
    796     std::string fn = mpsDir+"exmip1";
    797     m.readMps(fn.c_str(),"mps");
    798     ClpSimplex model;
    799     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    800                          m.getObjCoefficients(),
    801                          m.getRowLower(),m.getRowUpper());
    802     model.primal();
    803     int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
    804     double costIncrease[13];
    805     int sequenceIncrease[13];
    806     double costDecrease[13];
    807     int sequenceDecrease[13];
    808     // ranging
    809     model.dualRanging(13,which,costIncrease,sequenceIncrease,
    810                       costDecrease,sequenceDecrease);
    811     int i;
    812     for ( i=0;i<13;i++)
    813       printf("%d increase %g %d, decrease %g %d\n",
    814              i,costIncrease[i],sequenceIncrease[i],
    815              costDecrease[i],sequenceDecrease[i]);
    816     assert (fabs(costDecrease[3])<1.0e-4);
    817     assert (fabs(costIncrease[7]-1.0)<1.0e-4);
    818     model.setOptimizationDirection(-1);
    819     {
    820       int j;
    821       double * obj = model.objective();
    822       int n=model.numberColumns();
    823       for (j=0;j<n;j++)
    824         obj[j] *= -1.0;
    825     }
    826     double costIncrease2[13];
    827     int sequenceIncrease2[13];
    828     double costDecrease2[13];
    829     int sequenceDecrease2[13];
    830     // ranging
    831     model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
    832                       costDecrease2,sequenceDecrease2);
    833     for (i=0;i<13;i++) {
    834       assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
    835       assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
    836       assert (sequenceIncrease[i]==sequenceDecrease2[i]);
    837       assert (sequenceDecrease[i]==sequenceIncrease2[i]);
    838     }
    839     // Now delete all rows and see what happens
    840     model.deleteRows(model.numberRows(),which);
    841     model.primal();
    842     // ranging
    843     if (!model.dualRanging(8,which,costIncrease,sequenceIncrease,
    844                            costDecrease,sequenceDecrease)) {
    845       for (i=0;i<8;i++) {
    846         printf("%d increase %g %d, decrease %g %d\n",
    847                i,costIncrease[i],sequenceIncrease[i],
    848                costDecrease[i],sequenceDecrease[i]);
    849       }
    850     }
    851   }
    852   // Test primal ranging
    853   {   
    854     CoinMpsIO m;
    855     std::string fn = mpsDir+"exmip1";
    856     m.readMps(fn.c_str(),"mps");
    857     ClpSimplex model;
    858     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    859                          m.getObjCoefficients(),
    860                          m.getRowLower(),m.getRowUpper());
    861     model.primal();
    862     int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
    863     double valueIncrease[13];
    864     int sequenceIncrease[13];
    865     double valueDecrease[13];
    866     int sequenceDecrease[13];
    867     // ranging
    868     model.primalRanging(13,which,valueIncrease,sequenceIncrease,
    869                       valueDecrease,sequenceDecrease);
    870     int i;
    871     for ( i=0;i<13;i++)
    872       printf("%d increase %g %d, decrease %g %d\n",
    873              i,valueIncrease[i],sequenceIncrease[i],
    874              valueDecrease[i],sequenceDecrease[i]);
    875     assert (fabs(valueDecrease[3]-0.642857)<1.0e-4);
    876     assert (fabs(valueDecrease[8]-2.95113)<1.0e-4);
    877 #if 0
    878     // out until I find optimization bug
    879     // Test parametrics
    880     ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
    881     double rhs[]={ 1.0,2.0,3.0,4.0,5.0};
    882     double endingTheta=1.0;
    883     model2->scaling(0);
    884     model2->setLogLevel(63);
    885     model2->parametrics(0.0,endingTheta,0.1,
    886                         NULL,NULL,rhs,rhs,NULL);
    887 #endif
    888   }
    889   // Test binv etc
    890   {   
    891     /*
    892        Wolsey : Page 130
    893        max 4x1 -  x2
    894        7x1 - 2x2    <= 14
    895        x2    <= 3
    896        2x1 - 2x2    <= 3
    897        x1 in Z+, x2 >= 0
    898 
    899        note slacks are -1 in Clp so signs may be different
    900     */
    901    
    902     int n_cols = 2;
    903     int n_rows = 3;
    904    
    905     double obj[2] = {-4.0, 1.0};
    906     double collb[2] = {0.0, 0.0};
    907     double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
    908     double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
    909     double rowub[3] = {14.0, 3.0, 3.0};
    910    
    911     int rowIndices[5] =  {0,     2,    0,    1,    2};
    912     int colIndices[5] =  {0,     0,    1,    1,    1};
    913     double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
    914     CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
    915 
    916     ClpSimplex model;
    917     model.loadProblem(M, collb, colub, obj, rowlb, rowub);
    918     model.dual(0,1); // keep factorization
    919    
    920     //check that the tableau matches wolsey (B-1 A)
    921     // slacks in second part of binvA
    922     double * binvA = (double*) malloc((n_cols+n_rows) * sizeof(double));
    923    
    924     printf("B-1 A by row\n");
    925     int i;
    926     for( i = 0; i < n_rows; i++){
    927       model.getBInvARow(i, binvA,binvA+n_cols);
    928       printf("row: %d -> ",i);
    929       for(int j=0; j < n_cols+n_rows; j++){
    930         printf("%g, ", binvA[j]);
    931       }
    932       printf("\n");
    933     }
    934     // See if can re-use factorization AND arrays
    935     model.primal(0,3+4); // keep factorization
    936     // And do by column
    937     printf("B-1 A by column\n");
    938     for( i = 0; i < n_rows+n_cols; i++){
    939       model.getBInvACol(i, binvA);
    940       printf("column: %d -> ",i);
    941       for(int j=0; j < n_rows; j++){
    942         printf("%g, ", binvA[j]);
    943       }
    944       printf("\n");
    945     }
    946     free(binvA);
    947     model.setColUpper(1,2.0);
    948     model.dual(0,2+4); // use factorization and arrays
    949     model.dual(0,2); // hopefully will not use factorization
    950     model.primal(0,3+4); // keep factorization
    951     // but say basis has changed
    952     model.setWhatsChanged(model.whatsChanged()&(~512));
    953     model.dual(0,2); // hopefully will not use factorization
    954   }
    955   // test steepest edge
    956   {   
    957     CoinMpsIO m;
    958     std::string fn = mpsDir+"finnis";
    959     int returnCode = m.readMps(fn.c_str(),"mps");
    960     if (returnCode) {
    961       // probable cause is that gz not there
    962       fprintf(stderr,"Unable to open finnis.mps in %s!\n", mpsDir.c_str());
    963       fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
    964       fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
    965       exit(999);
    966     }
    967     ClpModel model;
    968     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
    969                     m.getColUpper(),
    970                     m.getObjCoefficients(),
    971                     m.getRowLower(),m.getRowUpper());
    972     ClpSimplex solution(model);
    973 
    974     solution.scaling(1);
    975     solution.setDualBound(1.0e8);
    976     //solution.factorization()->maximumPivots(1);
    977     //solution.setLogLevel(3);
    978     solution.setDualTolerance(1.0e-7);
    979     // set objective sense,
    980     ClpDualRowSteepest steep;
    981     solution.setDualRowPivotAlgorithm(steep);
    982     solution.setDblParam(ClpObjOffset,m.objectiveOffset());
    983     solution.dual();
    984   }
    985   // test normal solution
    986   {   
    987     CoinMpsIO m;
    988     std::string fn = mpsDir+"afiro";
    989     m.readMps(fn.c_str(),"mps");
    990     ClpSimplex solution;
    991     ClpModel model;
    992     // do twice - without and with scaling
    993     int iPass;
    994     for (iPass=0;iPass<2;iPass++) {
    995       // explicit row objective for testing
    996       int nr = m.getNumRows();
    997       double * rowObj = new double[nr];
    998       CoinFillN(rowObj,nr,0.0);
    999       model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1000                       m.getObjCoefficients(),
    1001                       m.getRowLower(),m.getRowUpper(),rowObj);
    1002       delete [] rowObj;
    1003       solution = ClpSimplex(model);
    1004       if (iPass) {
    1005         solution.scaling();
    1006       }
    1007       solution.dual();
    1008       solution.dual();
    1009       // test optimal
    1010       assert (solution.status()==0);
    1011       int numberColumns = solution.numberColumns();
    1012       int numberRows = solution.numberRows();
    1013       CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
    1014       double * objective = solution.objective();
    1015       double objValue = colsol.dotProduct(objective);
    1016       CoinRelFltEq eq(1.0e-8);
    1017       assert(eq(objValue,-4.6475314286e+02));
    1018       // Test auxiliary model
    1019       //solution.scaling(0);
    1020       solution.auxiliaryModel(63-2); // bounds may change
    1021       solution.dual();
    1022       solution.primal();
    1023       solution.allSlackBasis();
    1024       solution.dual();
    1025       assert(eq(solution.objectiveValue(),-4.6475314286e+02));
    1026       solution.auxiliaryModel(-1);
    1027       solution.dual();
    1028       assert(eq(solution.objectiveValue(),-4.6475314286e+02));
    1029       double * lower = solution.columnLower();
    1030       double * upper = solution.columnUpper();
    1031       double * sol = solution.primalColumnSolution();
    1032       double * result = new double[numberColumns];
    1033       CoinFillN ( result, numberColumns,0.0);
    1034       solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
    1035       int iRow , iColumn;
    1036       // see if feasible and dual feasible
    1037       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1038         double value = sol[iColumn];
    1039         assert(value<upper[iColumn]+1.0e-8);
    1040         assert(value>lower[iColumn]-1.0e-8);
    1041         value = objective[iColumn]-result[iColumn];
    1042         assert (value>-1.0e-5);
    1043         if (sol[iColumn]>1.0e-5)
    1044           assert (value<1.0e-5);
    1045       }
    1046       delete [] result;
    1047       result = new double[numberRows];
    1048       CoinFillN ( result, numberRows,0.0);
    1049       solution.matrix()->times(colsol, result);
    1050       lower = solution.rowLower();
    1051       upper = solution.rowUpper();
    1052       sol = solution.primalRowSolution();
    1053       for (iRow=0;iRow<numberRows;iRow++) {
    1054         double value = result[iRow];
    1055         assert(eq(value,sol[iRow]));
    1056         assert(value<upper[iRow]+1.0e-8);
    1057         assert(value>lower[iRow]-1.0e-8);
    1058       }
    1059       delete [] result;
    1060       // test row objective
    1061       double * rowObjective = solution.rowObjective();
    1062       CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
    1063       CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
    1064       // this sets up all slack basis
    1065       solution.createStatus();
    1066       solution.dual();
    1067       CoinFillN(rowObjective,numberRows,0.0);
    1068       CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
    1069       solution.dual();
    1070     }
    1071   }
    1072   // test unbounded
    1073   {   
    1074     CoinMpsIO m;
    1075     std::string fn = mpsDir+"brandy";
    1076     m.readMps(fn.c_str(),"mps");
    1077     ClpSimplex solution;
    1078     // do twice - without and with scaling
    1079     int iPass;
    1080     for (iPass=0;iPass<2;iPass++) {
    1081       solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1082                       m.getObjCoefficients(),
    1083                       m.getRowLower(),m.getRowUpper());
    1084       if (iPass)
    1085         solution.scaling();
    1086       solution.setOptimizationDirection(-1);
    1087       // test unbounded and ray
    1088 #ifdef DUAL
    1089       solution.setDualBound(100.0);
    1090       solution.dual();
    1091 #else
    1092       solution.primal();
    1093 #endif
    1094       assert (solution.status()==2);
    1095       int numberColumns = solution.numberColumns();
    1096       int numberRows = solution.numberRows();
    1097       double * lower = solution.columnLower();
    1098       double * upper = solution.columnUpper();
    1099       double * sol = solution.primalColumnSolution();
    1100       double * ray = solution.unboundedRay();
    1101       double * objective = solution.objective();
    1102       double objChange=0.0;
    1103       int iRow , iColumn;
    1104       // make sure feasible and columns form ray
    1105       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1106         double value = sol[iColumn];
    1107         assert(value<upper[iColumn]+1.0e-8);
    1108         assert(value>lower[iColumn]-1.0e-8);
    1109         value = ray[iColumn];
    1110         if (value>0.0)
    1111           assert(upper[iColumn]>1.0e30);
    1112         else if (value<0.0)
    1113           assert(lower[iColumn]<-1.0e30);
    1114         objChange += value*objective[iColumn];
    1115       }
    1116       // make sure increasing objective
    1117       assert(objChange>0.0);
    1118       double * result = new double[numberRows];
    1119       CoinFillN ( result, numberRows,0.0);
    1120       solution.matrix()->times(sol, result);
    1121       lower = solution.rowLower();
    1122       upper = solution.rowUpper();
    1123       sol = solution.primalRowSolution();
    1124       for (iRow=0;iRow<numberRows;iRow++) {
    1125         double value = result[iRow];
    1126         assert(eq(value,sol[iRow]));
    1127         assert(value<upper[iRow]+2.0e-8);
    1128         assert(value>lower[iRow]-2.0e-8);
    1129       }
    1130       CoinFillN ( result, numberRows,0.0);
    1131       solution.matrix()->times(ray, result);
    1132       // there may be small differences (especially if scaled)
    1133       for (iRow=0;iRow<numberRows;iRow++) {
    1134         double value = result[iRow];
    1135         if (value>1.0e-8)
    1136           assert(upper[iRow]>1.0e30);
    1137         else if (value<-1.0e-8)
    1138           assert(lower[iRow]<-1.0e30);
    1139       }
    1140       delete [] result;
    1141       delete [] ray;
    1142     }
    1143   }
    1144   // test infeasible
    1145   {   
    1146     CoinMpsIO m;
    1147     std::string fn = mpsDir+"brandy";
    1148     m.readMps(fn.c_str(),"mps");
    1149     ClpSimplex solution;
    1150     // do twice - without and with scaling
    1151     int iPass;
    1152     for (iPass=0;iPass<2;iPass++) {
    1153       solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1154                       m.getObjCoefficients(),
    1155                       m.getRowLower(),m.getRowUpper());
    1156       if (iPass)
    1157         solution.scaling();
    1158       // test infeasible and ray
    1159       solution.columnUpper()[0]=0.0;
    1160 #ifdef DUAL
    1161       solution.setDualBound(100.0);
    1162       solution.dual();
    1163 #else
    1164       solution.primal();
    1165 #endif
    1166       assert (solution.status()==1);
    1167       int numberColumns = solution.numberColumns();
    1168       int numberRows = solution.numberRows();
    1169       double * lower = solution.rowLower();
    1170       double * upper = solution.rowUpper();
    1171       double * ray = solution.infeasibilityRay();
    1172       assert(ray);
    1173       // construct proof of infeasibility
    1174       int iRow , iColumn;
    1175       double lo=0.0,up=0.0;
    1176       int nl=0,nu=0;
    1177       for (iRow=0;iRow<numberRows;iRow++) {
    1178         if (lower[iRow]>-1.0e20) {
    1179           lo += ray[iRow]*lower[iRow];
    1180         } else {
    1181           if (ray[iRow]>1.0e-8)
    1182             nl++;
    1183         }
    1184         if (upper[iRow]<1.0e20) {
    1185           up += ray[iRow]*upper[iRow];
    1186         } else {
    1187           if (ray[iRow]>1.0e-8)
    1188             nu++;
    1189         }
    1190       }
    1191       if (nl)
    1192         lo=-1.0e100;
    1193       if (nu)
    1194         up=1.0e100;
    1195       double * result = new double[numberColumns];
    1196       double lo2=0.0,up2=0.0;
    1197       CoinFillN ( result, numberColumns,0.0);
    1198       solution.matrix()->transposeTimes(ray, result);
    1199       lower = solution.columnLower();
    1200       upper = solution.columnUpper();
    1201       nl=nu=0;
    1202       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1203         if (result[iColumn]>1.0e-8) {
    1204           if (lower[iColumn]>-1.0e20)
    1205             lo2 += result[iColumn]*lower[iColumn];
    1206           else
    1207             nl++;
    1208           if (upper[iColumn]<1.0e20)
    1209             up2 += result[iColumn]*upper[iColumn];
    1210           else
    1211             nu++;
    1212         } else if (result[iColumn]<-1.0e-8) {
    1213           if (lower[iColumn]>-1.0e20)
    1214             up2 += result[iColumn]*lower[iColumn];
    1215           else
    1216             nu++;
    1217           if (upper[iColumn]<1.0e20)
    1218             lo2 += result[iColumn]*upper[iColumn];
    1219           else
    1220             nl++;
    1221         }
    1222       }
    1223       if (nl)
    1224         lo2=-1.0e100;
    1225       if (nu)
    1226         up2=1.0e100;
    1227       // make sure inconsistency
    1228       assert(lo2>up||up2<lo);
    1229       delete [] result;
    1230       delete [] ray;
    1231     }
    1232   }
    1233   // test delete and add
    1234   {   
    1235     CoinMpsIO m;
    1236     std::string fn = mpsDir+"brandy";
    1237     m.readMps(fn.c_str(),"mps");
    1238     ClpSimplex solution;
    1239     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1240                          m.getObjCoefficients(),
    1241                          m.getRowLower(),m.getRowUpper());
    1242     solution.dual();
    1243     CoinRelFltEq eq(1.0e-8);
    1244     assert(eq(solution.objectiveValue(),1.5185098965e+03));
    1245 
    1246     int numberColumns = solution.numberColumns();
    1247     int numberRows = solution.numberRows();
    1248     double * saveObj = new double [numberColumns];
    1249     double * saveLower = new double[numberRows+numberColumns];
    1250     double * saveUpper = new double[numberRows+numberColumns];
    1251     int * which = new int [numberRows+numberColumns];
    1252 
    1253     int numberElements = m.getMatrixByCol()->getNumElements();
    1254     int * starts = new int[numberRows+numberColumns];
    1255     int * index = new int[numberElements];
    1256     double * element = new double[numberElements];
    1257 
    1258     const CoinBigIndex * startM;
    1259     const int * lengthM;
    1260     const int * indexM;
    1261     const double * elementM;
    1262 
    1263     int n,nel;
    1264 
    1265     // delete non basic columns
    1266     n=0;
    1267     nel=0;
    1268     int iRow , iColumn;
    1269     const double * lower = m.getColLower();
    1270     const double * upper = m.getColUpper();
    1271     const double * objective = m.getObjCoefficients();
    1272     startM = m.getMatrixByCol()->getVectorStarts();
    1273     lengthM = m.getMatrixByCol()->getVectorLengths();
    1274     indexM = m.getMatrixByCol()->getIndices();
    1275     elementM = m.getMatrixByCol()->getElements();
    1276     starts[0]=0;
    1277     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1278       if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
    1279         saveObj[n]=objective[iColumn];
    1280         saveLower[n]=lower[iColumn];
    1281         saveUpper[n]=upper[iColumn];
    1282         int j;
    1283         for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
    1284           index[nel]=indexM[j];
    1285           element[nel++]=elementM[j];
    1286         }
    1287         which[n++]=iColumn;
    1288         starts[n]=nel;
    1289       }
    1290     }
    1291     solution.deleteColumns(n,which);
    1292     solution.dual();
    1293     // Put back
    1294     solution.addColumns(n,saveLower,saveUpper,saveObj,
    1295                         starts,index,element);
    1296     solution.dual();
    1297     assert(eq(solution.objectiveValue(),1.5185098965e+03));
    1298     // Delete all columns and add back
    1299     n=0;
    1300     nel=0;
    1301     starts[0]=0;
    1302     lower = m.getColLower();
    1303     upper = m.getColUpper();
    1304     objective = m.getObjCoefficients();
    1305     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1306       saveObj[n]=objective[iColumn];
    1307       saveLower[n]=lower[iColumn];
    1308       saveUpper[n]=upper[iColumn];
    1309       int j;
    1310       for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
    1311         index[nel]=indexM[j];
    1312         element[nel++]=elementM[j];
    1313       }
    1314       which[n++]=iColumn;
    1315       starts[n]=nel;
    1316     }
    1317     solution.deleteColumns(n,which);
    1318     solution.dual();
    1319     // Put back
    1320     solution.addColumns(n,saveLower,saveUpper,saveObj,
    1321                         starts,index,element);
    1322     solution.dual();
    1323     assert(eq(solution.objectiveValue(),1.5185098965e+03));
    1324 
    1325     // reload with original
    1326     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1327                          m.getObjCoefficients(),
    1328                          m.getRowLower(),m.getRowUpper());
    1329     // delete half rows
    1330     n=0;
    1331     nel=0;
    1332     lower = m.getRowLower();
    1333     upper = m.getRowUpper();
    1334     startM = m.getMatrixByRow()->getVectorStarts();
    1335     lengthM = m.getMatrixByRow()->getVectorLengths();
    1336     indexM = m.getMatrixByRow()->getIndices();
    1337     elementM = m.getMatrixByRow()->getElements();
    1338     starts[0]=0;
    1339     for (iRow=0;iRow<numberRows;iRow++) {
    1340       if ((iRow&1)==0) {
    1341         saveLower[n]=lower[iRow];
    1342         saveUpper[n]=upper[iRow];
    1343         int j;
    1344         for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
    1345           index[nel]=indexM[j];
    1346           element[nel++]=elementM[j];
    1347         }
    1348         which[n++]=iRow;
    1349         starts[n]=nel;
    1350       }
    1351     }
    1352     solution.deleteRows(n,which);
    1353     solution.dual();
    1354     // Put back
    1355     solution.addRows(n,saveLower,saveUpper,
    1356                         starts,index,element);
    1357     solution.dual();
    1358     assert(eq(solution.objectiveValue(),1.5185098965e+03));
    1359     solution.writeMps("yy.mps");
    1360     // Delete all rows
    1361     n=0;
    1362     nel=0;
    1363     lower = m.getRowLower();
    1364     upper = m.getRowUpper();
    1365     starts[0]=0;
    1366     for (iRow=0;iRow<numberRows;iRow++) {
    1367       saveLower[n]=lower[iRow];
    1368       saveUpper[n]=upper[iRow];
    1369       int j;
    1370       for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
    1371         index[nel]=indexM[j];
    1372         element[nel++]=elementM[j];
    1373       }
    1374       which[n++]=iRow;
    1375       starts[n]=nel;
    1376     }
    1377     solution.deleteRows(n,which);
    1378     solution.dual();
    1379     // Put back
    1380     solution.addRows(n,saveLower,saveUpper,
    1381                         starts,index,element);
    1382     solution.dual();
    1383     solution.writeMps("xx.mps");
    1384     assert(eq(solution.objectiveValue(),1.5185098965e+03));
    1385     // Zero out status array to give some interest
    1386     memset(solution.statusArray()+numberColumns,0,numberRows);
    1387     solution.primal(1);
    1388     assert(eq(solution.objectiveValue(),1.5185098965e+03));
    1389     // Delete all columns and rows
    1390     n=0;
    1391     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1392       which[n++]=iColumn;
    1393       starts[n]=nel;
    1394     }
    1395     solution.deleteColumns(n,which);
    1396     n=0;
    1397     for (iRow=0;iRow<numberRows;iRow++) {
    1398       which[n++]=iRow;
    1399       starts[n]=nel;
    1400     }
    1401     solution.deleteRows(n,which);
    1402 
    1403     delete [] saveObj;
    1404     delete [] saveLower;
    1405     delete [] saveUpper;
    1406     delete [] which;
    1407     delete [] starts;
    1408     delete [] index;
    1409     delete [] element;
    1410   }
    1411 #if 1
    1412   // Test barrier
    1413   {
    1414     CoinMpsIO m;
    1415     std::string fn = mpsDir+"exmip1";
    1416     m.readMps(fn.c_str(),"mps");
    1417     ClpInterior solution;
    1418     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1419                          m.getObjCoefficients(),
    1420                          m.getRowLower(),m.getRowUpper());
    1421     solution.primalDual();
    1422   }
    1423 #endif
    1424   // test network
    1425 #define QUADRATIC
    1426   if (1) {   
    1427     std::string fn = mpsDir+"input.130";
    1428     int numberColumns;
    1429     int numberRows;
    1430    
    1431     FILE * fp = fopen(fn.c_str(),"r");
    1432     if (!fp) {
    1433       // Try in Data
    1434       fn = "Data/Sample/input.130";
    1435       fp = fopen(fn.c_str(),"r");
    1436     }
    1437     if (!fp) {
    1438       fprintf(stderr,"Unable to open file input.130 in mpsDir or Data/Sample directory\n");
    1439     } else {
    1440       int problem;
    1441       char temp[100];
    1442       // read and skip
    1443       fscanf(fp,"%s",temp);
    1444       assert (!strcmp(temp,"BEGIN"));
    1445       fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows,
    1446              &numberColumns);
    1447       // scan down to SUPPLY
    1448       while (fgets(temp,100,fp)) {
    1449         if (!strncmp(temp,"SUPPLY",6))
    1450           break;
    1451       }
    1452       if (strncmp(temp,"SUPPLY",6)) {
    1453         fprintf(stderr,"Unable to find SUPPLY\n");
    1454         exit(2);
    1455       }
    1456       // get space for rhs
    1457       double * lower = new double[numberRows];
    1458       double * upper = new double[numberRows];
    1459       int i;
    1460       for (i=0;i<numberRows;i++) {
    1461         lower[i]=0.0;
    1462         upper[i]=0.0;
    1463       }
    1464       // ***** Remember to convert to C notation
    1465       while (fgets(temp,100,fp)) {
    1466         int row;
    1467         int value;
    1468         if (!strncmp(temp,"ARCS",4))
    1469           break;
    1470         sscanf(temp,"%d %d",&row,&value);
    1471         upper[row-1]=-value;
    1472         lower[row-1]=-value;
    1473       }
    1474       if (strncmp(temp,"ARCS",4)) {
    1475         fprintf(stderr,"Unable to find ARCS\n");
    1476         exit(2);
    1477       }
    1478       // number of columns may be underestimate
    1479       int * head = new int[2*numberColumns];
    1480       int * tail = new int[2*numberColumns];
    1481       int * ub = new int[2*numberColumns];
    1482       int * cost = new int[2*numberColumns];
    1483       // ***** Remember to convert to C notation
    1484       numberColumns=0;
    1485       while (fgets(temp,100,fp)) {
    1486         int iHead;
    1487         int iTail;
    1488         int iUb;
    1489         int iCost;
    1490         if (!strncmp(temp,"DEMAND",6))
    1491           break;
    1492         sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
    1493         iHead--;
    1494         iTail--;
    1495         head[numberColumns]=iHead;
    1496         tail[numberColumns]=iTail;
    1497         ub[numberColumns]=iUb;
    1498         cost[numberColumns]=iCost;
    1499         numberColumns++;
    1500       }
    1501       if (strncmp(temp,"DEMAND",6)) {
    1502         fprintf(stderr,"Unable to find DEMAND\n");
    1503         exit(2);
    1504       }
    1505       // ***** Remember to convert to C notation
    1506       while (fgets(temp,100,fp)) {
    1507         int row;
    1508         int value;
    1509         if (!strncmp(temp,"END",3))
    1510           break;
    1511         sscanf(temp,"%d %d",&row,&value);
    1512         upper[row-1]=value;
    1513         lower[row-1]=value;
    1514       }
    1515       if (strncmp(temp,"END",3)) {
    1516         fprintf(stderr,"Unable to find END\n");
    1517         exit(2);
    1518       }
    1519       printf("Problem %d has %d rows and %d columns\n",problem,
    1520              numberRows,numberColumns);
    1521       fclose(fp);
    1522       ClpSimplex  model;
    1523       // now build model
    1524      
    1525       double * objective =new double[numberColumns];
    1526       double * lowerColumn = new double[numberColumns];
    1527       double * upperColumn = new double[numberColumns];
    1528      
    1529       double * element = new double [2*numberColumns];
    1530       CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
    1531       int * row = new int[2*numberColumns];
    1532       start[numberColumns]=2*numberColumns;
    1533       for (i=0;i<numberColumns;i++) {
    1534         start[i]=2*i;
    1535         element[2*i]=-1.0;
    1536         element[2*i+1]=1.0;
    1537         row[2*i]=head[i];
    1538         row[2*i+1]=tail[i];
    1539         lowerColumn[i]=0.0;
    1540         upperColumn[i]=ub[i];
    1541         objective[i]=cost[i];
    1542       }
    1543       // Create Packed Matrix
    1544       CoinPackedMatrix matrix;
    1545       int * lengths = NULL;
    1546       matrix.assignMatrix(true,numberRows,numberColumns,
    1547                           2*numberColumns,element,row,start,lengths);
    1548       // load model
    1549       model.loadProblem(matrix,
    1550                         lowerColumn,upperColumn,objective,
    1551                         lower,upper);
    1552       model.factorization()->maximumPivots(200+model.numberRows()/100);
    1553       model.createStatus();
    1554       double time1 = CoinCpuTime();
    1555       model.dual();
    1556       std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    1557       ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix);
    1558       assert (plusMinus->getIndices()); // would be zero if not +- one
    1559       //ClpPlusMinusOneMatrix *plusminus_matrix;
    1560 
    1561       //plusminus_matrix = new ClpPlusMinusOneMatrix;
    1562 
    1563       //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
    1564       //                         plusMinus->startPositive(),plusMinus->startNegative());
    1565       model.loadProblem(*plusMinus,
    1566                         lowerColumn,upperColumn,objective,
    1567                         lower,upper);
    1568       //model.replaceMatrix( plusminus_matrix , true);
    1569       delete plusMinus;
    1570       //model.createStatus();
    1571       //model.initialSolve();
    1572       //model.writeMps("xx.mps");
    1573      
    1574       model.factorization()->maximumPivots(200+model.numberRows()/100);
    1575       model.createStatus();
    1576       time1 = CoinCpuTime();
    1577       model.dual();
    1578       std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    1579       ClpNetworkMatrix network(numberColumns,head,tail);
    1580       model.loadProblem(network,
    1581                         lowerColumn,upperColumn,objective,
    1582                         lower,upper);
    1583      
    1584       model.factorization()->maximumPivots(200+model.numberRows()/100);
    1585       model.createStatus();
    1586       time1 = CoinCpuTime();
    1587       model.dual();
    1588       std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    1589       delete [] lower;
    1590       delete [] upper;
    1591       delete [] head;
    1592       delete [] tail;
    1593       delete [] ub;
    1594       delete [] cost;
    1595       delete [] objective;
    1596       delete [] lowerColumn;
    1597       delete [] upperColumn;
    1598     }
    1599   }
    1600 #ifdef QUADRATIC
    1601   // Test quadratic to solve linear
    1602   if (1) {   
    1603     CoinMpsIO m;
    1604     std::string fn = mpsDir+"exmip1";
    1605     m.readMps(fn.c_str(),"mps");
    1606     ClpSimplex solution;
    1607     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1608                          m.getObjCoefficients(),
    1609                          m.getRowLower(),m.getRowUpper());
    1610     //solution.dual();
    1611     // get quadratic part
    1612     int numberColumns=solution.numberColumns();
    1613     int * start=new int [numberColumns+1];
    1614     int * column = new int[numberColumns];
    1615     double * element = new double[numberColumns];
    1616     int i;
    1617     start[0]=0;
    1618     int n=0;
    1619     int kk=numberColumns-1;
    1620     int kk2=numberColumns-1;
    1621     for (i=0;i<numberColumns;i++) {
    1622       if (i>=kk) {
    1623         column[n]=i;
    1624         if (i>=kk2)
    1625           element[n]=1.0e-1;
    1626         else
    1627           element[n]=0.0;
    1628         n++;
    1629       }
    1630       start[i+1]=n;
    1631     }
    1632     // Load up objective
    1633     solution.loadQuadraticObjective(numberColumns,start,column,element);
    1634     delete [] start;
    1635     delete [] column;
    1636     delete [] element;
    1637     //solution.quadraticSLP(50,1.0e-4);
    1638     double objValue = solution.getObjValue();
    1639     CoinRelFltEq eq(1.0e-4);
    1640     //assert(eq(objValue,820.0));
    1641     //solution.setLogLevel(63);
    1642     solution.primal();
    1643     int numberRows = solution.numberRows();
    1644 
    1645     double * rowPrimal = solution.primalRowSolution();
    1646     double * rowDual = solution.dualRowSolution();
    1647     double * rowLower = solution.rowLower();
    1648     double * rowUpper = solution.rowUpper();
    1649    
    1650     int iRow;
    1651     printf("Rows\n");
    1652     for (iRow=0;iRow<numberRows;iRow++) {
    1653       printf("%d primal %g dual %g low %g up %g\n",
    1654              iRow,rowPrimal[iRow],rowDual[iRow],
    1655              rowLower[iRow],rowUpper[iRow]);
    1656     }
    1657     double * columnPrimal = solution.primalColumnSolution();
    1658     double * columnDual = solution.dualColumnSolution();
    1659     double * columnLower = solution.columnLower();
    1660     double * columnUpper = solution.columnUpper();
    1661     objValue = solution.getObjValue();
    1662     int iColumn;
    1663     printf("Columns\n");
    1664     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1665       printf("%d primal %g dual %g low %g up %g\n",
    1666              iColumn,columnPrimal[iColumn],columnDual[iColumn],
    1667              columnLower[iColumn],columnUpper[iColumn]);
    1668     }
    1669     //assert(eq(objValue,3.2368421));
    1670     //exit(77);
    1671   }
    1672   // Test quadratic
    1673   if (1) {   
    1674     CoinMpsIO m;
    1675     std::string fn = mpsDir+"share2qp";
    1676     //fn = "share2qpb";
    1677     m.readMps(fn.c_str(),"mps");
    1678     ClpSimplex model;
    1679     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1680                          m.getObjCoefficients(),
    1681                          m.getRowLower(),m.getRowUpper());
    1682     model.dual();
    1683     // get quadratic part
    1684     int * start=NULL;
    1685     int * column = NULL;
    1686     double * element = NULL;
    1687     m.readQuadraticMps(NULL,start,column,element,2);
    1688     int column2[200];
    1689     double element2[200];
    1690     int start2[80];
    1691     int j;
    1692     start2[0]=0;
    1693     int nel=0;
    1694     bool good=false;
    1695     for (j=0;j<79;j++) {
    1696       if (start[j]==start[j+1]) {
    1697         column2[nel]=j;
    1698         element2[nel]=0.0;
    1699         nel++;
    1700       } else {
    1701         int i;
    1702         for (i=start[j];i<start[j+1];i++) {
    1703           column2[nel]=column[i];
    1704           element2[nel++]=element[i];
    1705         }
    1706       }
    1707       start2[j+1]=nel;
    1708     }
    1709     // Load up objective
    1710     if (good)
    1711       model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
    1712     else
    1713       model.loadQuadraticObjective(model.numberColumns(),start,column,element);
    1714     delete [] start;
    1715     delete [] column;
    1716     delete [] element;
    1717     int numberColumns=model.numberColumns();
    1718 #if 0
    1719     model.nonlinearSLP(50,1.0e-4);
    1720 #else
    1721     // Get feasible
    1722     ClpObjective * saveObjective = model.objectiveAsObject()->clone();
    1723     ClpLinearObjective zeroObjective(NULL,numberColumns);
    1724     model.setObjective(&zeroObjective);
    1725     model.dual();
    1726     model.setObjective(saveObjective);
    1727     delete saveObjective;
    1728 #endif
    1729     //model.setLogLevel(63);
    1730     //exit(77);
    1731     model.setFactorizationFrequency(10);
    1732     model.primal();
    1733     double objValue = model.getObjValue();
    1734     CoinRelFltEq eq(1.0e-4);
    1735     assert(eq(objValue,-400.92));
    1736   }
    1737   if (0) {   
    1738     CoinMpsIO m;
    1739     std::string fn = "./beale";
    1740     //fn = "./jensen";
    1741     m.readMps(fn.c_str(),"mps");
    1742     ClpSimplex solution;
    1743     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
    1744                          m.getObjCoefficients(),
    1745                          m.getRowLower(),m.getRowUpper());
    1746     solution.setDblParam(ClpObjOffset,m.objectiveOffset());
    1747     solution.dual();
    1748     // get quadratic part
    1749     int * start=NULL;
    1750     int * column = NULL;
    1751     double * element = NULL;
    1752     m.readQuadraticMps(NULL,start,column,element,2);
    1753     // Load up objective
    1754     solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
    1755     delete [] start;
    1756     delete [] column;
    1757     delete [] element;
    1758     solution.primal(1);
    1759     solution.nonlinearSLP(50,1.0e-4);
    1760     double objValue = solution.getObjValue();
    1761     CoinRelFltEq eq(1.0e-4);
    1762     assert(eq(objValue,0.5));
    1763     solution.primal();
    1764     objValue = solution.getObjValue();
    1765     assert(eq(objValue,0.5));
    1766   }
    1767 #endif 
    1768 }
    1769 void CbcClpUnitTest (const CbcModel & saveModel, std::string& miplibDir)
     29//#############################################################################
     30
     31void CbcClpUnitTest (const CbcModel & saveModel, std::string& miplibDir,
     32                     bool unitTestOnly)
    177033{
    177134  unsigned int m ;
     
    177336  FILE * fp;
    177437  bool doTest=false;
    1775 //  const char dirsep =  CoinFindDirSeparator();
    1776  
     38
    177739  // Set directory containing miplib data files.
    177840  std::string test1 = miplibDir +"p0033";
     
    181779  strategy.push_back(zz_strategy_zz) ; \
    181880  objValue.push_back(zz_objValue_zz) ;
    1819  
    1820   /*
    1821     Load up the problem vector. Note that the row counts here include the
    1822     objective function.
    1823    
    1824   */
    1825   // 0 for no test, 1 for some, 2 for many, 3 for all
     81
     82  if (unitTestOnly) {
     83    PUSH_MPS("p0033",16,33,3089,2520.57,7);
     84    PUSH_MPS("p0201",133,201,7615,6875.0,7);
     85    PUSH_MPS("p0548",176,548,8691,315.29,7);
     86  } else {
     87    /*
     88      Load up the problem vector. Note that the row counts here include the
     89      objective function.
     90    */
     91    // 0 for no test, 1 for some, 2 for many, 3 for all
    182692#define HOWMANY 2
    182793#if HOWMANY
    182894#if HOWMANY>1
    1829   PUSH_MPS("10teams",230,2025,924,917,7);
    1830 #endif
    1831   PUSH_MPS("air03",124,10757,340160,338864.25,7);
    1832 #if HOWMANY==3
    1833   PUSH_MPS("air04",823,8904,56137,55535.436,8);
    1834   PUSH_MPS("air05",426,7195,26374,25877.609,8);
     95    PUSH_MPS("10teams",230,2025,924,917,7);
     96#endif
     97    PUSH_MPS("air03",124,10757,340160,338864.25,7);
     98#if HOWMANY==3
     99    PUSH_MPS("air04",823,8904,56137,55535.436,8);
     100    PUSH_MPS("air05",426,7195,26374,25877.609,8);
    1835101#endif
    1836102  //    PUSH_MPS("arki001",1048,1388,7580813.0459,7579599.80787,7);
    1837   PUSH_MPS("bell3a",123,133,878430.32,862578.64,7);
    1838 #if HOWMANY>1
    1839   PUSH_MPS("bell5",91,104,8966406.49,8608417.95,7);
    1840 #endif
    1841   PUSH_MPS("blend2",274,353,7.598985,6.9156751140,7);
    1842 #if HOWMANY>1
    1843   PUSH_MPS("cap6000",2176,6000,-2451377,-2451537.325,7);
    1844 #endif
    1845   //    PUSH_MPS("dano3mip",3202,13873,728.1111,576.23162474,7);
    1846   //PUSH_MPS("danoint",664,521,65.67,62.637280418,7);
    1847   PUSH_MPS("dcmulti",290,548,188182,183975.5397,7);
    1848   PUSH_MPS("dsbmip",1182,1886,-305.19817501,-305.19817501,7);
    1849   PUSH_MPS("egout",98,141,568.101,149.589,7);
    1850   PUSH_MPS("enigma",21,100,0.0,0.0,7);
    1851 #if HOWMANY==3
    1852   PUSH_MPS("fast0507",507,63009,174,172.14556668,7);
    1853 #endif
    1854   PUSH_MPS("fiber",363,1298,405935.18000,156082.51759,7);
    1855 #if HOWMANY>1
    1856   PUSH_MPS("fixnet6",478,878,3983,1200.88,7);
    1857 #endif
    1858   PUSH_MPS("flugpl",18,18,1201500,1167185.7,7);
    1859   PUSH_MPS("gen",780,870,112313,112130.0,7);
    1860 #if HOWMANY>1
    1861   PUSH_MPS("gesa2",1392,1224,25779856.372,25476489.678,7);
    1862   PUSH_MPS("gesa2_o",1248,1224,25779856.372,25476489.678,7);
    1863 #endif
    1864   PUSH_MPS("gesa3",1368,1152,27991042.648,27833632.451,7);
    1865   PUSH_MPS("gesa3_o",1224,1152,27991042.648,27833632.451,7);
    1866   PUSH_MPS("gt2",29,188,21166.000,13460.233074,7);
    1867 #if HOWMANY==3
    1868   PUSH_MPS("harp2",112,2993,-73899798.00,-74353341.502,7);
    1869 #endif
    1870   PUSH_MPS("khb05250",101,1350,106940226,95919464.0,7);
    1871 #if HOWMANY>1
    1872   PUSH_MPS("l152lav",97,1989,4722,4656.36,7);
    1873 #endif
    1874   PUSH_MPS("lseu",28,89,1120,834.68,7);
    1875   PUSH_MPS("misc03",96,160,3360,1910.,7);
    1876   PUSH_MPS("misc06",820,1808,12850.8607,12841.6,7);
    1877 #if HOWMANY>1
    1878   PUSH_MPS("misc07",212,260,2810,1415.0,7);
    1879   PUSH_MPS("mitre",2054,10724,115155,114740.5184,7);
    1880 #endif
    1881   PUSH_MPS("mod008",6,319,307,290.9,7);
    1882   PUSH_MPS("mod010",146,2655,6548,6532.08,7);
    1883 #if HOWMANY==3
    1884   PUSH_MPS("mod011",4480,10958,-54558535,-62121982.55,7);
    1885   PUSH_MPS("modglob",291,422,20740508,20430947.,7);
    1886   PUSH_MPS("noswot",182,128,-43,-43.0,7);
    1887 #endif
    1888 #if HOWMANY>1
    1889   PUSH_MPS("nw04",36,87482,16862,16310.66667,7);
    1890 #endif
    1891   PUSH_MPS("p0033",16,33,3089,2520.57,7);
    1892   PUSH_MPS("p0201",133,201,7615,6875.0,7);
    1893   PUSH_MPS("p0282",241,282,258411,176867.50,7);
    1894   PUSH_MPS("p0548",176,548,8691,315.29,7);
    1895   PUSH_MPS("p2756",755,2756,3124,2688.75,7);
    1896 #if HOWMANY==3
    1897   PUSH_MPS("pk1",45,86,11.0,0.0,7);
    1898 #endif
    1899 #if HOWMANY>1
    1900   PUSH_MPS("pp08a",136,240,7350.0,2748.3452381,7);
    1901   PUSH_MPS("pp08aCUTS",246,240,7350.0,5480.6061563,7);
    1902 #endif
    1903 #if HOWMANY==3
    1904   PUSH_MPS("qiu",1192,840,-132.873137,-931.638857,7);
    1905 #endif
    1906   PUSH_MPS("qnet1",503,1541,16029.692681,14274.102667,7);
    1907   PUSH_MPS("qnet1_o",456,1541,16029.692681,12095.571667,7);
    1908   PUSH_MPS("rentacar",6803,9557,30356761,28806137.644,7);
    1909   PUSH_MPS("rgn",24,180,82.1999,48.7999,7);
    1910 #if HOWMANY==3
    1911   PUSH_MPS("rout",291,556,1077.56,981.86428571,7);
    1912   PUSH_MPS("set1ch",492,712,54537.75,32007.73,7);
    1913 #endif
    1914   //    PUSH_MPS("seymour",4944,1372,423,403.84647413,7);
    1915   PUSH_MPS("stein27",118,27,18,13.0,7);
    1916 #if HOWMANY>1
    1917   PUSH_MPS("stein45",331,45,30,22.0,7);
    1918 #endif
    1919   PUSH_MPS("vpm1",234,378,20,15.4167,7);
    1920   PUSH_MPS("vpm2",234,378,13.75,9.8892645972,7);
    1921 #endif
     103    PUSH_MPS("bell3a",123,133,878430.32,862578.64,7);
     104#if HOWMANY>1
     105    PUSH_MPS("bell5",91,104,8966406.49,8608417.95,7);
     106#endif
     107    PUSH_MPS("blend2",274,353,7.598985,6.9156751140,7);
     108#if HOWMANY>1
     109    PUSH_MPS("cap6000",2176,6000,-2451377,-2451537.325,7);
     110#endif
     111    //    PUSH_MPS("dano3mip",3202,13873,728.1111,576.23162474,7);
     112    //PUSH_MPS("danoint",664,521,65.67,62.637280418,7);
     113    PUSH_MPS("dcmulti",290,548,188182,183975.5397,7);
     114    PUSH_MPS("dsbmip",1182,1886,-305.19817501,-305.19817501,7);
     115    PUSH_MPS("egout",98,141,568.101,149.589,7);
     116    PUSH_MPS("enigma",21,100,0.0,0.0,7);
     117#if HOWMANY==3
     118    PUSH_MPS("fast0507",507,63009,174,172.14556668,7);
     119#endif
     120    PUSH_MPS("fiber",363,1298,405935.18000,156082.51759,7);
     121#if HOWMANY>1
     122    PUSH_MPS("fixnet6",478,878,3983,1200.88,7);
     123#endif
     124    PUSH_MPS("flugpl",18,18,1201500,1167185.7,7);
     125    PUSH_MPS("gen",780,870,112313,112130.0,7);
     126#if HOWMANY>1
     127    PUSH_MPS("gesa2",1392,1224,25779856.372,25476489.678,7);
     128    PUSH_MPS("gesa2_o",1248,1224,25779856.372,25476489.678,7);
     129#endif
     130    PUSH_MPS("gesa3",1368,1152,27991042.648,27833632.451,7);
     131    PUSH_MPS("gesa3_o",1224,1152,27991042.648,27833632.451,7);
     132    PUSH_MPS("gt2",29,188,21166.000,13460.233074,7);
     133#if HOWMANY==3
     134    PUSH_MPS("harp2",112,2993,-73899798.00,-74353341.502,7);
     135#endif
     136    PUSH_MPS("khb05250",101,1350,106940226,95919464.0,7);
     137#if HOWMANY>1
     138    PUSH_MPS("l152lav",97,1989,4722,4656.36,7);
     139#endif
     140    PUSH_MPS("lseu",28,89,1120,834.68,7);
     141    PUSH_MPS("misc03",96,160,3360,1910.,7);
     142    PUSH_MPS("misc06",820,1808,12850.8607,12841.6,7);
     143#if HOWMANY>1
     144    PUSH_MPS("misc07",212,260,2810,1415.0,7);
     145    PUSH_MPS("mitre",2054,10724,115155,114740.5184,7);
     146#endif
     147    PUSH_MPS("mod008",6,319,307,290.9,7);
     148    PUSH_MPS("mod010",146,2655,6548,6532.08,7);
     149#if HOWMANY==3
     150    PUSH_MPS("mod011",4480,10958,-54558535,-62121982.55,7);
     151    PUSH_MPS("modglob",291,422,20740508,20430947.,7);
     152    PUSH_MPS("noswot",182,128,-43,-43.0,7);
     153#endif
     154#if HOWMANY>1
     155    PUSH_MPS("nw04",36,87482,16862,16310.66667,7);
     156#endif
     157    PUSH_MPS("p0033",16,33,3089,2520.57,7);
     158    PUSH_MPS("p0201",133,201,7615,6875.0,7);
     159    PUSH_MPS("p0282",241,282,258411,176867.50,7);
     160    PUSH_MPS("p0548",176,548,8691,315.29,7);
     161    PUSH_MPS("p2756",755,2756,3124,2688.75,7);
     162#if HOWMANY==3
     163    PUSH_MPS("pk1",45,86,11.0,0.0,7);
     164#endif
     165#if HOWMANY>1
     166    PUSH_MPS("pp08a",136,240,7350.0,2748.3452381,7);
     167    PUSH_MPS("pp08aCUTS",246,240,7350.0,5480.6061563,7);
     168#endif
     169#if HOWMANY==3
     170    PUSH_MPS("qiu",1192,840,-132.873137,-931.638857,7);
     171#endif
     172    PUSH_MPS("qnet1",503,1541,16029.692681,14274.102667,7);
     173    PUSH_MPS("qnet1_o",456,1541,16029.692681,12095.571667,7);
     174    PUSH_MPS("rentacar",6803,9557,30356761,28806137.644,7);
     175    PUSH_MPS("rgn",24,180,82.1999,48.7999,7);
     176#if HOWMANY==3
     177    PUSH_MPS("rout",291,556,1077.56,981.86428571,7);
     178    PUSH_MPS("set1ch",492,712,54537.75,32007.73,7);
     179#endif
     180    //    PUSH_MPS("seymour",4944,1372,423,403.84647413,7);
     181    PUSH_MPS("stein27",118,27,18,13.0,7);
     182#if HOWMANY>1
     183    PUSH_MPS("stein45",331,45,30,22.0,7);
     184#endif
     185    PUSH_MPS("vpm1",234,378,20,15.4167,7);
     186    PUSH_MPS("vpm2",234,378,13.75,9.8892645972,7);
     187#endif
     188  }
    1922189#undef PUSH_MPS
    1923190   
     
    1929196  */
    1930197  for (m = 0 ; m < mpsName.size() ; m++) {
    1931     std::cerr << "  processing mps file: " << mpsName[m]
    1932               << " (" << m+1 << " out of " << mpsName.size() << ")" << std::endl ;
     198    std::cout << "  processing mps file: " << mpsName[m]
     199              << " (" << m+1 << " out of " << mpsName.size() << ")\n";
    1933200    /*
    1934201      Stage 1: Read the MPS
     
    1939206    std::string fn = miplibDir+mpsName[m] ;
    1940207    model->solver()->readMps(fn.c_str(),"") ;
    1941     int nr = model->getNumRows() ;
    1942     int nc = model->getNumCols() ;
    1943     assert(nr == nRows[m]) ;
    1944     assert(nc == nCols[m]) ;
     208    assert(model->getNumRows() == nRows[m]) ;
     209    assert(model->getNumCols() == nCols[m]) ;
     210
    1945211    /*
    1946212      Stage 2: Call solver to solve the problem.
    1947      
    1948213      then check the return code and objective.
    1949      
    1950214    */
    1951215
    1952216    double startTime = CoinCpuTime()+CoinCpuTimeJustChildren();
    1953     if (model->getMaximumNodes()>200000)
     217    if (model->getMaximumNodes()>200000) {
    1954218      model->setMaximumNodes(200000);
     219    }
    1955220    OsiClpSolverInterface * si =
    1956221      dynamic_cast<OsiClpSolverInterface *>(model->solver()) ;
     
    1979244          below *= multiplier;
    1980245        }
    1981         if (above<1.0e12)
     246        if (above<1.0e12) {
    1982247          largest = CoinMax(largest,above);
    1983         if (below<1.0e12)
     248        }
     249        if (below<1.0e12) {
    1984250          largest = CoinMax(largest,below);
     251        }
    1985252      }
    1986253     
     
    2000267          below *= multiplier;
    2001268        }
    2002         if (above<1.0e12)
     269        if (above<1.0e12) {
    2003270          largest = CoinMax(largest,above);
    2004         if (below<1.0e12)
     271        }
     272        if (below<1.0e12) {
    2005273          largest = CoinMax(largest,below);
     274        }
    2006275      }
    2007276      //std::cout<<"Largest (scaled) away from bound "<<largest<<std::endl;
    2008       modelC->setDualBound(CoinMax(1.0001e8,CoinMin(1000.0*largest,1.00001e10)));
     277      modelC->setDualBound(CoinMax(1.0001e8,
     278                                   CoinMin(1000.0*largest,1.00001e10)));
    2009279    }
    2010280    model->setMinimumDrop(min(5.0e-2,
    2011                                  fabs(model->getMinimizationObjValue())*1.0e-3+1.0e-4));
    2012     if (model->getNumCols()<500)
     281                              fabs(model->getMinimizationObjValue())*1.0e-3+1.0e-4));
     282    if (model->getNumCols()<500) {
    2013283      model->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
    2014     else if (model->getNumCols()<5000)
     284    } else if (model->getNumCols()<5000) {
    2015285      model->setMaximumCutPassesAtRoot(100); // use minimum drop
    2016     else
     286    } else {
    2017287      model->setMaximumCutPassesAtRoot(20);
     288    }
    2018289    // If defaults then increase trust for small models
    2019290    if (model->numberStrong()==5&&model->numberBeforeTrust()==10) {
    2020291      int numberColumns = model->getNumCols();
    2021       if (numberColumns<=50)
     292      if (numberColumns<=50) {
    2022293        model->setNumberBeforeTrust(1000);
    2023       else if (numberColumns<=100)
     294      } else if (numberColumns<=100) {
    2024295        model->setNumberBeforeTrust(100);
    2025       else if (numberColumns<=300)
     296      } else if (numberColumns<=300) {
    2026297        model->setNumberBeforeTrust(50);
     298      }
    2027299    }
    2028300    model->branchAndBound();
     
    2048320      CoinRelFltEq eq(1.0e-3) ;
    2049321      if (eq(soln,objValue[m])) {
    2050         std::cerr
     322        std::cout
    2051323          <<"cbc_clp"<<" "
    2052324          << soln << " = " << objValue[m] << " ; okay";
    2053325        numProbSolved++;
    2054326      } else  {
    2055         std::cerr <<"cbc_clp" <<" " <<soln << " != " <<objValue[m] << "; error=" ;
    2056         std::cerr <<fabs(objValue[m] - soln);
     327        std::cout <<"cbc_clp" <<" " <<soln << " != " <<objValue[m]
     328                  << "; error=" << fabs(objValue[m] - soln);
    2057329      }
    2058330    } else {
    2059       std::cerr << "error; too many nodes" ;
    2060     }
    2061     std::cerr<<" - took " <<timeOfSolution<<" seconds."<<std::endl;
     331      std::cout << "error; too many nodes" ;
     332    }
     333    std::cout<<" - took " <<timeOfSolution<<" seconds.\n";
    2062334    timeTaken += timeOfSolution;
    2063335    delete model;
    2064336  }
    2065   std::cerr
     337  std::cout
    2066338    <<"cbc_clp"
    2067339    <<" solved "
Note: See TracChangeset for help on using the changeset viewer.