source: trunk/Cbc/src/CoinSolve.cpp @ 312

Last change on this file since 312 was 312, checked in by andreasw, 14 years ago

finished (for now) conversion of Cbc to autotools

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 131.2 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3   
4#include "CoinPragma.hpp"
5
6#include <cassert>
7#include <cstdio>
8#include <cmath>
9#include <cfloat>
10#include <cstring>
11#include <iostream>
12
13
14#include "CoinPragma.hpp"
15#include "CoinHelperFunctions.hpp"
16// Same version as CBC
17#define CBCVERSION "1.01.00"
18
19#include "CoinMpsIO.hpp"
20
21#include "ClpFactorization.hpp"
22#include "CoinTime.hpp"
23#include "ClpSimplex.hpp"
24#include "ClpSimplexOther.hpp"
25#include "ClpSolve.hpp"
26#include "ClpPackedMatrix.hpp"
27#include "ClpPlusMinusOneMatrix.hpp"
28#include "ClpNetworkMatrix.hpp"
29#include "ClpDualRowSteepest.hpp"
30#include "ClpDualRowDantzig.hpp"
31#include "ClpLinearObjective.hpp"
32#include "ClpPrimalColumnSteepest.hpp"
33#include "ClpPrimalColumnDantzig.hpp"
34#include "ClpPresolve.hpp"
35#include "CbcOrClpParam.hpp"
36#include "OsiRowCutDebugger.hpp"
37#ifdef DMALLOC
38#include "dmalloc.h"
39#endif
40#ifdef WSSMP_BARRIER
41#define FOREIGN_BARRIER
42#endif
43#ifdef UFL_BARRIER
44#define FOREIGN_BARRIER
45#endif
46#ifdef TAUCS_BARRIER
47#define FOREIGN_BARRIER
48#endif
49#include "CoinWarmStartBasis.hpp"
50
51#include "OsiSolverInterface.hpp"
52#include "OsiCuts.hpp"
53#include "OsiRowCut.hpp"
54#include "OsiColCut.hpp"
55
56#include "CglPreProcess.hpp"
57#include "CglCutGenerator.hpp"
58#include "CglGomory.hpp"
59#include "CglProbing.hpp"
60#include "CglKnapsackCover.hpp"
61#include "CglRedSplit.hpp"
62#include "CglClique.hpp"
63#include "CglFlowCover.hpp"
64#include "CglMixedIntegerRounding2.hpp"
65#include "CglTwomir.hpp"
66
67#include "CbcModel.hpp"
68#include "CbcHeuristic.hpp"
69#include "CbcHeuristicLocal.hpp"
70#include "CbcHeuristicGreedy.hpp"
71#include "CbcHeuristicFPump.hpp"
72#include "CbcTreeLocal.hpp"
73#include "CbcCompareActual.hpp"
74#include "CbcBranchActual.hpp"
75#include  "CbcOrClpParam.hpp"
76#include  "CbcCutGenerator.hpp"
77#include  "CbcStrategy.hpp"
78
79#include "OsiClpSolverInterface.hpp"
80#ifdef CBC_AMPL
81#include "Cbc_ampl.h"
82static bool usingAmpl=false;
83#endif
84static double totalTime=0.0;
85static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
86static bool maskMatches(std::string & mask, std::string & check);
87//#############################################################################
88
89#ifdef NDEBUG
90#undef NDEBUG
91#endif
92// Allow for interrupts
93// But is this threadsafe ? (so switched off by option)
94
95#include "CoinSignal.hpp"
96static CbcModel * currentBranchModel = NULL;
97
98extern "C" {
99   static void signal_handler(int whichSignal)
100   {
101      if (currentBranchModel!=NULL) 
102         currentBranchModel->setMaximumNodes(0); // stop at next node
103      return;
104   }
105}
106
107int mainTest (int argc, const char *argv[],int algorithm,
108              ClpSimplex empty, bool doPresolve,int switchOff);
109void CbcClpUnitTest (const CbcModel & saveModel);
110int CbcOrClpRead_mode=1;
111FILE * CbcOrClpReadCommand=stdin;
112static bool noPrinting=false;
113static int * analyze(OsiClpSolverInterface * solverMod, int & numberChanged, double & increment,
114                     bool changeInt)
115{
116  OsiSolverInterface * solver = solverMod->clone();
117  if (0) {
118    // just get increment
119    CbcModel model(*solver);
120    model.analyzeObjective();
121    double increment2=model.getCutoffIncrement();
122    printf("initial cutoff increment %g\n",increment2);
123  }
124  const double *objective = solver->getObjCoefficients() ;
125  const double *lower = solver->getColLower() ;
126  const double *upper = solver->getColUpper() ;
127  int numberColumns = solver->getNumCols() ;
128  int numberRows = solver->getNumRows();
129  double direction = solver->getObjSense();
130  int iRow,iColumn;
131
132  // Row copy
133  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
134  const double * elementByRow = matrixByRow.getElements();
135  const int * column = matrixByRow.getIndices();
136  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
137  const int * rowLength = matrixByRow.getVectorLengths();
138
139  // Column copy
140  CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
141  const double * element = matrixByCol.getElements();
142  const int * row = matrixByCol.getIndices();
143  const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
144  const int * columnLength = matrixByCol.getVectorLengths();
145
146  const double * rowLower = solver->getRowLower();
147  const double * rowUpper = solver->getRowUpper();
148
149  char * ignore = new char [numberRows];
150  int * changed = new int[numberColumns];
151  int * which = new int[numberRows];
152  double * changeRhs = new double[numberRows];
153  memset(changeRhs,0,numberRows*sizeof(double));
154  memset(ignore,0,numberRows);
155  numberChanged=0;
156  int numberInteger=0;
157  for (iColumn=0;iColumn<numberColumns;iColumn++) {
158    if (upper[iColumn] > lower[iColumn]+1.0e-8&&solver->isInteger(iColumn)) 
159      numberInteger++;
160  }
161  bool finished=false;
162  while (!finished) {
163    int saveNumberChanged = numberChanged;
164    for (iRow=0;iRow<numberRows;iRow++) {
165      int numberContinuous=0;
166      double value1=0.0,value2=0.0;
167      bool allIntegerCoeff=true;
168      double sumFixed=0.0;
169      int jColumn1=-1,jColumn2=-1;
170      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
171        int jColumn = column[j];
172        double value = elementByRow[j];
173        if (upper[jColumn] > lower[jColumn]+1.0e-8) {
174          if (!solver->isInteger(jColumn)) {
175            if (numberContinuous==0) {
176              jColumn1=jColumn;
177              value1=value;
178            } else {
179              jColumn2=jColumn;
180              value2=value;
181            }
182            numberContinuous++;
183          } else {
184            if (fabs(value-floor(value+0.5))>1.0e-12)
185              allIntegerCoeff=false;
186          }
187        } else {
188          sumFixed += lower[jColumn]*value;
189        }
190      }
191      double low = rowLower[iRow];
192      if (low>-1.0e20) {
193        low -= sumFixed;
194        if (fabs(low-floor(low+0.5))>1.0e-12)
195          allIntegerCoeff=false;
196      }
197      double up = rowUpper[iRow];
198      if (up<1.0e20) {
199        up -= sumFixed;
200        if (fabs(up-floor(up+0.5))>1.0e-12)
201          allIntegerCoeff=false;
202      }
203      if (!allIntegerCoeff)
204        continue; // can't do
205      if (numberContinuous==1) {
206        // see if really integer
207        // This does not allow for complicated cases
208        if (low==up) {
209          if (fabs(value1)>1.0e-3) {
210            value1 = 1.0/value1;
211            if (fabs(value1-floor(value1+0.5))<1.0e-12) {
212              // integer
213              changed[numberChanged++]=jColumn1;
214              solver->setInteger(jColumn1);
215              if (upper[jColumn1]>1.0e20)
216                solver->setColUpper(jColumn1,1.0e20);
217              if (lower[jColumn1]<-1.0e20)
218                solver->setColLower(jColumn1,-1.0e20);
219            }
220          }
221        } else {
222          if (fabs(value1)>1.0e-3) {
223            value1 = 1.0/value1;
224            if (fabs(value1-floor(value1+0.5))<1.0e-12) {
225              // This constraint will not stop it being integer
226              ignore[iRow]=1;
227            }
228          }
229        }
230      } else if (numberContinuous==2) {
231        if (low==up) {
232          /* need general theory - for now just look at 2 cases -
233             1 - +- 1 one in column and just costs i.e. matching objective
234             2 - +- 1 two in column but feeds into G/L row which will try and minimize
235          */
236          if (fabs(value1)==1.0&&value1*value2==-1.0&&!lower[jColumn1]
237              &&!lower[jColumn2]) {
238            int n=0;
239            int i;
240            double objChange=direction*(objective[jColumn1]+objective[jColumn2]);
241            double bound = CoinMin(upper[jColumn1],upper[jColumn2]);
242            bound = CoinMin(bound,1.0e20);
243            for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
244              int jRow = row[i];
245              double value = element[i];
246              if (jRow!=iRow) {
247                which[n++]=jRow;
248                changeRhs[jRow]=value;
249              }
250            }
251            for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
252              int jRow = row[i];
253              double value = element[i];
254              if (jRow!=iRow) {
255                if (!changeRhs[jRow]) {
256                  which[n++]=jRow;
257                  changeRhs[jRow]=value;
258                } else {
259                  changeRhs[jRow]+=value;
260                }
261              }
262            }
263            if (objChange>=0.0) {
264              // see if all rows OK
265              bool good=true;
266              for (i=0;i<n;i++) {
267                int jRow = which[i];
268                double value = changeRhs[jRow];
269                if (value) {
270                  value *= bound;
271                  if (rowLength[jRow]==1) {
272                    if (value>0.0) {
273                      double rhs = rowLower[jRow];
274                      if (rhs>0.0) {
275                        double ratio =rhs/value;
276                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
277                          good=false;
278                      }
279                    } else {
280                      double rhs = rowUpper[jRow];
281                      if (rhs<0.0) {
282                        double ratio =rhs/value;
283                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
284                          good=false;
285                      }
286                    }
287                  } else if (rowLength[jRow]==2) {
288                    if (value>0.0) {
289                      if (rowLower[jRow]>-1.0e20)
290                        good=false;
291                    } else {
292                      if (rowUpper[jRow]<1.0e20)
293                        good=false;
294                    }
295                  } else {
296                    good=false;
297                  }
298                }
299              }
300              if (good) {
301                // both can be integer
302                changed[numberChanged++]=jColumn1;
303                solver->setInteger(jColumn1);
304                if (upper[jColumn1]>1.0e20)
305                  solver->setColUpper(jColumn1,1.0e20);
306                if (lower[jColumn1]<-1.0e20)
307                  solver->setColLower(jColumn1,-1.0e20);
308                changed[numberChanged++]=jColumn2;
309                solver->setInteger(jColumn2);
310                if (upper[jColumn2]>1.0e20)
311                  solver->setColUpper(jColumn2,1.0e20);
312                if (lower[jColumn2]<-1.0e20)
313                  solver->setColLower(jColumn2,-1.0e20);
314              }
315            }
316            // clear
317            for (i=0;i<n;i++) {
318              changeRhs[which[i]]=0.0;
319            }
320          }
321        }
322      }
323    }
324    for (iColumn=0;iColumn<numberColumns;iColumn++) {
325      if (upper[iColumn] > lower[iColumn]+1.0e-8&&!solver->isInteger(iColumn)) {
326        double value;
327        value = upper[iColumn];
328        if (value<1.0e20&&fabs(value-floor(value+0.5))>1.0e-12) 
329          continue;
330        value = lower[iColumn];
331        if (value>-1.0e20&&fabs(value-floor(value+0.5))>1.0e-12) 
332          continue;
333        bool integer=true;
334        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
335          int iRow = row[j];
336          if (!ignore[iRow]) {
337            integer=false;
338            break;
339          }
340        }
341        if (integer) {
342          // integer
343          changed[numberChanged++]=iColumn;
344          solver->setInteger(iColumn);
345          if (upper[iColumn]>1.0e20)
346            solver->setColUpper(iColumn,1.0e20);
347          if (lower[iColumn]<-1.0e20)
348            solver->setColLower(iColumn,-1.0e20);
349        }
350      }
351    }
352    finished = numberChanged==saveNumberChanged;
353  }
354  delete [] which;
355  delete [] changeRhs;
356  delete [] ignore;
357  if (numberInteger&&!noPrinting)
358    printf("%d integer variables",numberInteger);
359  if (changeInt) {
360    if (!noPrinting) {
361      if (numberChanged)
362        printf(" and %d variables made integer\n",numberChanged);
363      else
364        printf("\n");
365    }
366    delete [] ignore;
367    //increment=0.0;
368    if (!numberChanged) {
369      delete [] changed;
370      delete solver;
371      return NULL;
372    } else {
373      for (iColumn=0;iColumn<numberColumns;iColumn++) {
374        if (solver->isInteger(iColumn))
375          solverMod->setInteger(iColumn);
376      }
377      delete solver;
378      return changed;
379    }
380  } else {
381    if (!noPrinting) {
382      if (numberChanged)
383        printf(" and %d variables could be made integer\n",numberChanged);
384      else
385        printf("\n");
386    }
387    // just get increment
388    CbcModel model(*solver);
389    if (noPrinting)
390      model.setLogLevel(0);
391    model.analyzeObjective();
392    double increment2=model.getCutoffIncrement();
393    if (increment2>increment) {
394      if (!noPrinting)
395        printf("cutoff increment increased from %g to %g\n",increment,increment2);
396      increment=increment2;
397    }
398    delete solver;
399    numberChanged=0;
400    delete [] changed;
401    return NULL;
402  }
403}
404int main (int argc, const char *argv[])
405{
406  /* Note
407     This is meant as a stand-alone executable to do as much of coin as possible.
408     It should only have one solver known to it.
409  */
410  {
411    double time1 = CoinCpuTime(),time2;
412    bool goodModel=false;
413    CoinSighandler_t saveSignal=SIG_DFL;
414    // register signal handler
415    saveSignal = signal(SIGINT,signal_handler);
416    // Set up all non-standard stuff
417    OsiClpSolverInterface solver1;
418    CbcModel model(solver1);
419    CbcModel * babModel = NULL;
420    model.setNumberBeforeTrust(21);
421    int cutPass=-1234567;
422    OsiSolverInterface * solver = model.solver();
423    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
424    ClpSimplex * lpSolver = clpSolver->getModelPtr();
425    clpSolver->messageHandler()->setLogLevel(0) ;
426    model.messageHandler()->setLogLevel(1);
427    // For priorities etc
428    int * priorities=NULL;
429    int * branchDirection=NULL;
430    double * pseudoDown=NULL;
431    double * pseudoUp=NULL;
432#ifdef CBC_AMPL
433    ampl_info info;
434    if (argc>2&&!strcmp(argv[2],"-AMPL")) {
435      usingAmpl=true;
436      int returnCode = readAmpl(&info,argc,const_cast<char **>(argv));
437      if (returnCode)
438        return returnCode;
439      CbcOrClpRead_mode=2; // so will start with parameters
440      // see if log in list
441      noPrinting=true;
442      for (int i=1;i<info.numberArguments;i++) {
443        if (!strcmp(info.arguments[i],"log")) {
444          if (i<info.numberArguments-1&&atoi(info.arguments[i+1])>0)
445            noPrinting=false;
446          break;
447        }
448      }
449      if (noPrinting) {
450        model.messageHandler()->setLogLevel(0);
451        setCbcOrClpPrinting(false);
452      }
453      if (!noPrinting)
454        printf("%d rows, %d columns and %d elements\n",
455               info.numberRows,info.numberColumns,info.numberElements);
456      solver->loadProblem(info.numberColumns,info.numberRows,info.starts,
457                          info.rows,info.elements,
458                          info.columnLower,info.columnUpper,info.objective,
459                          info.rowLower,info.rowUpper);
460      // If we had a solution use it
461      if (info.primalSolution) {
462        solver->setColSolution(info.primalSolution);
463      }
464      // status
465      if (info.rowStatus) {
466        unsigned char * statusArray = lpSolver->statusArray();
467        int i;
468        for (i=0;i<info.numberColumns;i++)
469          statusArray[i]=(char)info.columnStatus[i];
470        statusArray+=info.numberColumns;
471        for (i=0;i<info.numberRows;i++)
472          statusArray[i]=(char)info.rowStatus[i];
473        CoinWarmStartBasis * basis = lpSolver->getBasis();
474        solver->setWarmStart(basis);
475        delete basis;
476      }
477      freeArrays1(&info);
478      // modify objective if necessary
479      solver->setObjSense(info.direction);
480      solver->setDblParam(OsiObjOffset,info.offset);
481      // Set integer variables
482      for (int i=info.numberColumns-info.numberBinary-info.numberIntegers;
483           i<info.numberColumns;i++)
484        solver->setInteger(i);
485      goodModel=true;
486      // change argc etc
487      argc = info.numberArguments;
488      argv = const_cast<const char **>(info.arguments);
489    }
490#endif   
491    // default action on import
492    int allowImportErrors=0;
493    int keepImportNames=1;
494    int doIdiot=-1;
495    int outputFormat=2;
496    int slpValue=-1;
497    int printOptions=0;
498    int printMode=0;
499    int presolveOptions=0;
500    int substitution=3;
501    int dualize=0;
502    int doCrash=0;
503    int doSprint=-1;
504    int doScaling=1;
505    // set reasonable defaults
506    int preSolve=5;
507    int preProcess=4;
508    bool preSolveFile=false;
509   
510    double djFix=1.0e100;
511    double gapRatio=1.0e100;
512    double tightenFactor=0.0;
513    lpSolver->setPerturbation(50);
514    lpSolver->messageHandler()->setPrefix(false);
515    const char dirsep =  CoinFindDirSeparator();
516    std::string directory = (dirsep == '/' ? "./" : ".\\");
517    std::string defaultDirectory = directory;
518    std::string importFile ="";
519    std::string exportFile ="default.mps";
520    std::string importBasisFile ="";
521    std::string importPriorityFile ="";
522    std::string debugFile="";
523    std::string printMask="";
524    double * debugValues = NULL;
525    int numberDebugValues = -1;
526    int basisHasValues=0;
527    std::string exportBasisFile ="default.bas";
528    std::string saveFile ="default.prob";
529    std::string restoreFile ="default.prob";
530    std::string solutionFile ="stdout";
531    std::string solutionSaveFile ="solution.file";
532#define CBCMAXPARAMETERS 200
533    CbcOrClpParam parameters[CBCMAXPARAMETERS];
534    int numberParameters ;
535    establishParams(numberParameters,parameters) ;
536    parameters[whichParam(BASISIN,numberParameters,parameters)].setStringValue(importBasisFile);
537    parameters[whichParam(PRIORITYIN,numberParameters,parameters)].setStringValue(importPriorityFile);
538    parameters[whichParam(BASISOUT,numberParameters,parameters)].setStringValue(exportBasisFile);
539    parameters[whichParam(DEBUG,numberParameters,parameters)].setStringValue(debugFile);
540    parameters[whichParam(PRINTMASK,numberParameters,parameters)].setStringValue(printMask);
541    parameters[whichParam(DIRECTORY,numberParameters,parameters)].setStringValue(directory);
542    parameters[whichParam(DUALBOUND,numberParameters,parameters)].setDoubleValue(lpSolver->dualBound());
543    parameters[whichParam(DUALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->dualTolerance());
544    parameters[whichParam(EXPORT,numberParameters,parameters)].setStringValue(exportFile);
545    parameters[whichParam(IDIOT,numberParameters,parameters)].setIntValue(doIdiot);
546    parameters[whichParam(IMPORT,numberParameters,parameters)].setStringValue(importFile);
547    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].setDoubleValue(1.0e-8);
548    int slog = whichParam(SOLVERLOGLEVEL,numberParameters,parameters);
549    int log = whichParam(LOGLEVEL,numberParameters,parameters);
550    parameters[slog].setIntValue(0);
551    parameters[log].setIntValue(1);
552    parameters[whichParam(MAXFACTOR,numberParameters,parameters)].setIntValue(lpSolver->factorizationFrequency());
553    parameters[whichParam(MAXITERATION,numberParameters,parameters)].setIntValue(lpSolver->maximumIterations());
554    parameters[whichParam(OUTPUTFORMAT,numberParameters,parameters)].setIntValue(outputFormat);
555    parameters[whichParam(PRESOLVEPASS,numberParameters,parameters)].setIntValue(preSolve);
556    parameters[whichParam(PERTVALUE,numberParameters,parameters)].setIntValue(lpSolver->perturbation());
557    parameters[whichParam(PRIMALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->primalTolerance());
558    parameters[whichParam(PRIMALWEIGHT,numberParameters,parameters)].setDoubleValue(lpSolver->infeasibilityCost());
559    parameters[whichParam(RESTORE,numberParameters,parameters)].setStringValue(restoreFile);
560    parameters[whichParam(SAVE,numberParameters,parameters)].setStringValue(saveFile);
561    //parameters[whichParam(TIMELIMIT,numberParameters,parameters)].setDoubleValue(1.0e8);
562    parameters[whichParam(TIMELIMIT_BAB,numberParameters,parameters)].setDoubleValue(1.0e8);
563    parameters[whichParam(SOLUTION,numberParameters,parameters)].setStringValue(solutionFile);
564    parameters[whichParam(SAVESOL,numberParameters,parameters)].setStringValue(solutionSaveFile);
565    parameters[whichParam(SPRINT,numberParameters,parameters)].setIntValue(doSprint);
566    parameters[whichParam(SUBSTITUTION,numberParameters,parameters)].setIntValue(substitution);
567    parameters[whichParam(DUALIZE,numberParameters,parameters)].setIntValue(dualize);
568    model.setNumberBeforeTrust(5);
569    parameters[whichParam(NUMBERBEFORE,numberParameters,parameters)].setIntValue(5);
570    parameters[whichParam(MAXNODES,numberParameters,parameters)].setIntValue(model.getMaximumNodes());
571    model.setNumberStrong(5);
572    parameters[whichParam(STRONGBRANCHING,numberParameters,parameters)].setIntValue(model.numberStrong());
573    parameters[whichParam(INFEASIBILITYWEIGHT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcInfeasibilityWeight));
574    parameters[whichParam(INTEGERTOLERANCE,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcIntegerTolerance));
575    parameters[whichParam(INCREMENT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcCutoffIncrement));
576    // Set up likely cut generators and defaults
577    parameters[whichParam(PREPROCESS,numberParameters,parameters)].setCurrentOption("sos");
578    parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].setIntValue(128|64|1);
579    parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].setIntValue(-1);
580    parameters[whichParam(MAXHOTITS,numberParameters,parameters)].setIntValue(100);
581    parameters[whichParam(CUTSSTRATEGY,numberParameters,parameters)].setCurrentOption("on");
582    parameters[whichParam(HEURISTICSTRATEGY,numberParameters,parameters)].setCurrentOption("on");
583    int verbose=0;
584    CglGomory gomoryGen;
585    // try larger limit
586    gomoryGen.setLimitAtRoot(512);
587    gomoryGen.setLimit(50);
588    // set default action (0=off,1=on,2=root)
589    int gomoryAction=3;
590    parameters[whichParam(GOMORYCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
591
592    CglProbing probingGen;
593    probingGen.setUsingObjective(true);
594    probingGen.setMaxPass(3);
595    probingGen.setMaxPassRoot(3);
596    // Number of unsatisfied variables to look at
597    probingGen.setMaxProbe(10);
598    probingGen.setMaxProbeRoot(50);
599    // How far to follow the consequences
600    probingGen.setMaxLook(10);
601    probingGen.setMaxLookRoot(50);
602    probingGen.setMaxLookRoot(10);
603    // Only look at rows with fewer than this number of elements
604    probingGen.setMaxElements(200);
605    probingGen.setRowCuts(3);
606    // set default action (0=off,1=on,2=root)
607    int probingAction=3;
608    parameters[whichParam(PROBINGCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
609
610    CglKnapsackCover knapsackGen;
611    //knapsackGen.switchOnExpensive();
612    // set default action (0=off,1=on,2=root)
613    int knapsackAction=3;
614    parameters[whichParam(KNAPSACKCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
615
616    CglRedSplit redsplitGen;
617    //redsplitGen.setLimit(100);
618    // set default action (0=off,1=on,2=root)
619    // Off as seems to give some bad cuts
620    int redsplitAction=2;
621    parameters[whichParam(REDSPLITCUTS,numberParameters,parameters)].setCurrentOption("root");
622
623    CglClique cliqueGen(false,true);
624    cliqueGen.setStarCliqueReport(false);
625    cliqueGen.setRowCliqueReport(false);
626    cliqueGen.setMinViolation(0.1);
627    // set default action (0=off,1=on,2=root)
628    int cliqueAction=3;
629    parameters[whichParam(CLIQUECUTS,numberParameters,parameters)].setCurrentOption("ifmove");
630
631    CglMixedIntegerRounding2 mixedGen;
632    // set default action (0=off,1=on,2=root)
633    int mixedAction=3;
634    parameters[whichParam(MIXEDCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
635
636    CglFlowCover flowGen;
637    // set default action (0=off,1=on,2=root)
638    int flowAction=3;
639    parameters[whichParam(FLOWCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
640
641    CglTwomir twomirGen;
642    twomirGen.setMaxElements(250);
643    // set default action (0=off,1=on,2=root)
644    int twomirAction=2;
645    parameters[whichParam(TWOMIRCUTS,numberParameters,parameters)].setCurrentOption("root");
646
647    bool useRounding=true;
648    parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption("on");
649    bool useFpump=true;
650    parameters[whichParam(FPUMP,numberParameters,parameters)].setCurrentOption("on");
651    bool useGreedy=true;
652    parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption("on");
653    bool useCombine=true;
654    parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption("on");
655    bool useLocalTree=false;
656    parameters[whichParam(COSTSTRATEGY,numberParameters,parameters)].setCurrentOption("off");
657    bool useCosts=false;
658   
659    // total number of commands read
660    int numberGoodCommands=0;
661    // Set false if user does anything advanced
662    bool defaultSettings=true;
663
664    // Hidden stuff for barrier
665    int choleskyType = 0;
666    int gamma=0;
667    int scaleBarrier=0;
668    int doKKT=0;
669    int crossover=2; // do crossover unless quadratic
670    // For names
671    int lengthName = 0;
672    std::vector<std::string> rowNames;
673    std::vector<std::string> columnNames;
674   
675    std::string field;
676    if (!noPrinting) {
677      std::cout<<"Coin Cbc and Clp Solver version "<<CBCVERSION
678               <<", build "<<__DATE__<<std::endl;
679      // Print command line
680      if (argc>1) {
681        printf("command line - ");
682        for (int i=0;i<argc;i++)
683          printf("%s ",argv[i]);
684        printf("\n");
685      }
686    }
687    while (1) {
688      // next command
689      field=CoinReadGetCommand(argc,argv);
690      // exit if null or similar
691      if (!field.length()) {
692        if (numberGoodCommands==1&&goodModel) {
693          // we just had file name - do branch and bound
694          field="branch";
695        } else if (!numberGoodCommands) {
696          // let's give the sucker a hint
697          std::cout
698            <<"CoinSolver takes input from arguments ( - switches to stdin)"
699            <<std::endl
700            <<"Enter ? for list of commands or help"<<std::endl;
701          field="-";
702        } else {
703          break;
704        }
705      }
706     
707      // see if ? at end
708      int numberQuery=0;
709      if (field!="?"&&field!="???") {
710        int length = field.length();
711        int i;
712        for (i=length-1;i>0;i--) {
713          if (field[i]=='?') 
714            numberQuery++;
715          else
716            break;
717        }
718        field=field.substr(0,length-numberQuery);
719      }
720      // find out if valid command
721      int iParam;
722      int numberMatches=0;
723      int firstMatch=-1;
724      for ( iParam=0; iParam<numberParameters; iParam++ ) {
725        int match = parameters[iParam].matches(field);
726        if (match==1) {
727          numberMatches = 1;
728          firstMatch=iParam;
729          break;
730        } else {
731          if (match&&firstMatch<0)
732            firstMatch=iParam;
733          numberMatches += match>>1;
734        }
735      }
736      if (iParam<numberParameters&&!numberQuery) {
737        // found
738        CbcOrClpParam found = parameters[iParam];
739        CbcOrClpParameterType type = found.type();
740        int valid;
741        numberGoodCommands++;
742        if (type==BAB&&goodModel) {
743          // check if any integers
744          if (!lpSolver->integerInformation())
745            type=DUALSIMPLEX;
746        }
747        if (type==GENERALQUERY) {
748#ifdef CBC_AMPL
749          if (verbose<4&&usingAmpl)
750            verbose +=4;
751#endif
752          if (verbose<4) {
753            std::cout<<"In argument list keywords have leading - "
754              ", -stdin or just - switches to stdin"<<std::endl;
755            std::cout<<"One command per line (and no -)"<<std::endl;
756            std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
757            std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
758            std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
759            std::cout<<"abcd value sets value"<<std::endl;
760            std::cout<<"Commands are:"<<std::endl;
761          } else {
762            std::cout<<"Cbc options are set within AMPL with commands like:"<<std::endl<<std::endl;
763            std::cout<<"         option cbc_options \"cuts=root log=2 feas=on slog=1\""<<std::endl<<std::endl;
764            std::cout<<"only maximize, dual, primal, help and quit are recognized without ="<<std::endl;
765          }
766          int maxAcross=5;
767          if ((verbose%4)!=0)
768            maxAcross=1;
769          int limits[]={1,51,101,151,201,251,301,351,401};
770          std::vector<std::string> types;
771          types.push_back("Double parameters:");
772          types.push_back("Branch and Cut double parameters:");
773          types.push_back("Integer parameters:");
774          types.push_back("Branch and Cut integer parameters:");
775          types.push_back("Keyword parameters:");
776          types.push_back("Branch and Cut keyword parameters:");
777          types.push_back("Actions or string parameters:");
778          types.push_back("Branch and Cut actions:");
779          int iType;
780          for (iType=0;iType<8;iType++) {
781            int across=0;
782            if ((verbose%4)!=0)
783              std::cout<<std::endl;
784            std::cout<<types[iType]<<std::endl;
785            if ((verbose&2)!=0)
786              std::cout<<std::endl;
787            for ( iParam=0; iParam<numberParameters; iParam++ ) {
788              int type = parameters[iParam].type();
789              if (parameters[iParam].displayThis()&&type>=limits[iType]
790                  &&type<limits[iType+1]) {
791                // but skip if not useful for ampl (and in ampl mode)
792                if (verbose>=4&&(parameters[iParam].whereUsed()&4)==0)
793                  continue;
794                if (!across) {
795                  if ((verbose&2)==0) 
796                    std::cout<<"  ";
797                  else
798                    std::cout<<"Command ";
799                }
800                std::cout<<parameters[iParam].matchName()<<"  ";
801                across++;
802                if (across==maxAcross) {
803                  across=0;
804                  if ((verbose%4)!=0) {
805                    // put out description as well
806                    if ((verbose&1)!=0) 
807                      std::cout<<parameters[iParam].shortHelp();
808                    std::cout<<std::endl;
809                    if ((verbose&2)!=0) {
810                      std::cout<<"---- description"<<std::endl;
811                      parameters[iParam].printLongHelp();
812                      std::cout<<"----"<<std::endl<<std::endl;
813                    }
814                  } else {
815                    std::cout<<std::endl;
816                  }
817                }
818              }
819            }
820            if (across)
821              std::cout<<std::endl;
822          }
823        } else if (type==FULLGENERALQUERY) {
824          std::cout<<"Full list of commands is:"<<std::endl;
825          int maxAcross=5;
826          int limits[]={1,51,101,151,201,251,301,351,401};
827          std::vector<std::string> types;
828          types.push_back("Double parameters:");
829          types.push_back("Branch and Cut double parameters:");
830          types.push_back("Integer parameters:");
831          types.push_back("Branch and Cut integer parameters:");
832          types.push_back("Keyword parameters:");
833          types.push_back("Branch and Cut keyword parameters:");
834          types.push_back("Actions or string parameters:");
835          types.push_back("Branch and Cut actions:");
836          int iType;
837          for (iType=0;iType<8;iType++) {
838            int across=0;
839            std::cout<<types[iType]<<"  ";
840            for ( iParam=0; iParam<numberParameters; iParam++ ) {
841              int type = parameters[iParam].type();
842              if (type>=limits[iType]
843                  &&type<limits[iType+1]) {
844                if (!across)
845                  std::cout<<"  ";
846                std::cout<<parameters[iParam].matchName()<<"  ";
847                across++;
848                if (across==maxAcross) {
849                  std::cout<<std::endl;
850                  across=0;
851                }
852              }
853            }
854            if (across)
855              std::cout<<std::endl;
856          }
857        } else if (type<101) {
858          // get next field as double
859          double value = CoinReadGetDoubleField(argc,argv,&valid);
860          if (!valid) {
861            if (type<51) {
862              parameters[iParam].setDoubleParameter(lpSolver,value);
863            } else if (type<81) {
864              parameters[iParam].setDoubleParameter(model,value);
865            } else {
866              parameters[iParam].setDoubleParameter(lpSolver,value);
867              switch(type) {
868              case DJFIX:
869                djFix=value;
870                preSolve=5;
871                defaultSettings=false; // user knows what she is doing
872                break;
873              case GAPRATIO:
874                gapRatio=value;
875                break;
876              case TIGHTENFACTOR:
877                tightenFactor=value;
878                defaultSettings=false; // user knows what she is doing
879                break;
880              default:
881                abort();
882              }
883            }
884          } else if (valid==1) {
885            abort();
886          } else {
887            std::cout<<parameters[iParam].name()<<" has value "<<
888              parameters[iParam].doubleValue()<<std::endl;
889          }
890        } else if (type<201) {
891          // get next field as int
892          int value = CoinReadGetIntField(argc,argv,&valid);
893          if (!valid) {
894            if (type<151) {
895              if (parameters[iParam].type()==PRESOLVEPASS)
896                preSolve = value;
897              else if (parameters[iParam].type()==IDIOT)
898                doIdiot = value;
899              else if (parameters[iParam].type()==SPRINT)
900                doSprint = value;
901              else if (parameters[iParam].type()==OUTPUTFORMAT)
902                outputFormat = value;
903              else if (parameters[iParam].type()==SLPVALUE)
904                slpValue = value;
905              else if (parameters[iParam].type()==PRESOLVEOPTIONS)
906                presolveOptions = value;
907              else if (parameters[iParam].type()==PRINTOPTIONS)
908                printOptions = value;
909              else if (parameters[iParam].type()==SUBSTITUTION)
910                substitution = value;
911              else if (parameters[iParam].type()==DUALIZE)
912                dualize = value;
913              else if (parameters[iParam].type()==CUTPASS)
914                cutPass = value;
915              else if (parameters[iParam].type()==VERBOSE)
916                verbose = value;
917              else if (parameters[iParam].type()==FPUMPITS)
918                { useFpump = true;parameters[iParam].setIntValue(value);}
919              parameters[iParam].setIntParameter(lpSolver,value);
920            } else {
921              parameters[iParam].setIntParameter(model,value);
922            }
923          } else if (valid==1) {
924            abort();
925          } else {
926            std::cout<<parameters[iParam].name()<<" has value "<<
927              parameters[iParam].intValue()<<std::endl;
928          }
929        } else if (type<301) {
930          // one of several strings
931          std::string value = CoinReadGetString(argc,argv);
932          int action = parameters[iParam].parameterOption(value);
933          if (action<0) {
934            if (value!="EOL") {
935              // no match
936              parameters[iParam].printOptions();
937            } else {
938              // print current value
939              std::cout<<parameters[iParam].name()<<" has value "<<
940                parameters[iParam].currentOption()<<std::endl;
941            }
942          } else {
943            parameters[iParam].setCurrentOption(action,!noPrinting);
944            // for now hard wired
945            switch (type) {
946            case DIRECTION:
947              if (action==0)
948                lpSolver->setOptimizationDirection(1);
949              else if (action==1)
950                lpSolver->setOptimizationDirection(-1);
951              else
952                lpSolver->setOptimizationDirection(0);
953              break;
954            case DUALPIVOT:
955              if (action==0) {
956                ClpDualRowSteepest steep(3);
957                lpSolver->setDualRowPivotAlgorithm(steep);
958              } else if (action==1) {
959                //ClpDualRowDantzig dantzig;
960                ClpDualRowSteepest dantzig(5);
961                lpSolver->setDualRowPivotAlgorithm(dantzig);
962              } else if (action==2) {
963                // partial steep
964                ClpDualRowSteepest steep(2);
965                lpSolver->setDualRowPivotAlgorithm(steep);
966              } else {
967                ClpDualRowSteepest steep;
968                lpSolver->setDualRowPivotAlgorithm(steep);
969              }
970              break;
971            case PRIMALPIVOT:
972              if (action==0) {
973                ClpPrimalColumnSteepest steep(3);
974                lpSolver->setPrimalColumnPivotAlgorithm(steep);
975              } else if (action==1) {
976                ClpPrimalColumnSteepest steep(0);
977                lpSolver->setPrimalColumnPivotAlgorithm(steep);
978              } else if (action==2) {
979                ClpPrimalColumnDantzig dantzig;
980                lpSolver->setPrimalColumnPivotAlgorithm(dantzig);
981              } else if (action==3) {
982                ClpPrimalColumnSteepest steep(2);
983                lpSolver->setPrimalColumnPivotAlgorithm(steep);
984              } else if (action==4) {
985                ClpPrimalColumnSteepest steep(1);
986                lpSolver->setPrimalColumnPivotAlgorithm(steep);
987              } else if (action==5) {
988                ClpPrimalColumnSteepest steep(4);
989                lpSolver->setPrimalColumnPivotAlgorithm(steep);
990              } else if (action==6) {
991                ClpPrimalColumnSteepest steep(10);
992                lpSolver->setPrimalColumnPivotAlgorithm(steep);
993              }
994              break;
995            case SCALING:
996              lpSolver->scaling(action);
997              solver->setHintParam(OsiDoScale,action!=0,OsiHintTry);
998              doScaling = 1-action;
999              break;
1000            case AUTOSCALE:
1001              lpSolver->setAutomaticScaling(action!=0);
1002              break;
1003            case SPARSEFACTOR:
1004              lpSolver->setSparseFactorization((1-action)!=0);
1005              break;
1006            case BIASLU:
1007              lpSolver->factorization()->setBiasLU(action);
1008              break;
1009            case PERTURBATION:
1010              if (action==0)
1011                lpSolver->setPerturbation(50);
1012              else
1013                lpSolver->setPerturbation(100);
1014              break;
1015            case ERRORSALLOWED:
1016              allowImportErrors = action;
1017              break;
1018            case INTPRINT:
1019              printMode=action;
1020              break;
1021              //case ALGORITHM:
1022              //algorithm  = action;
1023              //defaultSettings=false; // user knows what she is doing
1024              //abort();
1025              //break;
1026            case KEEPNAMES:
1027              keepImportNames = 1-action;
1028              break;
1029            case PRESOLVE:
1030              if (action==0)
1031                preSolve = 5;
1032              else if (action==1)
1033                preSolve=0;
1034              else if (action==2)
1035                preSolve=10;
1036              else
1037                preSolveFile=true;
1038              break;
1039            case PFI:
1040              lpSolver->factorization()->setForrestTomlin(action==0);
1041              break;
1042            case CRASH:
1043              doCrash=action;
1044              break;
1045            case MESSAGES:
1046              lpSolver->messageHandler()->setPrefix(action!=0);
1047              break;
1048            case CHOLESKY:
1049              choleskyType = action;
1050              break;
1051            case GAMMA:
1052              gamma=action;
1053              break;
1054            case BARRIERSCALE:
1055              scaleBarrier=action;
1056              break;
1057            case KKT:
1058              doKKT=action;
1059              break;
1060            case CROSSOVER:
1061              crossover=action;
1062              break;
1063            case GOMORYCUTS:
1064              defaultSettings=false; // user knows what she is doing
1065              gomoryAction = action;
1066              break;
1067            case PROBINGCUTS:
1068              defaultSettings=false; // user knows what she is doing
1069              probingAction = action;
1070              break;
1071            case KNAPSACKCUTS:
1072              defaultSettings=false; // user knows what she is doing
1073              knapsackAction = action;
1074              break;
1075            case REDSPLITCUTS:
1076              defaultSettings=false; // user knows what she is doing
1077              redsplitAction = action;
1078              break;
1079            case CLIQUECUTS:
1080              defaultSettings=false; // user knows what she is doing
1081              cliqueAction = action;
1082              break;
1083            case FLOWCUTS:
1084              defaultSettings=false; // user knows what she is doing
1085              flowAction = action;
1086              break;
1087            case MIXEDCUTS:
1088              defaultSettings=false; // user knows what she is doing
1089              mixedAction = action;
1090              break;
1091            case TWOMIRCUTS:
1092              defaultSettings=false; // user knows what she is doing
1093              twomirAction = action;
1094              break;
1095            case ROUNDING:
1096              defaultSettings=false; // user knows what she is doing
1097              useRounding = action;
1098              break;
1099            case FPUMP:
1100              defaultSettings=false; // user knows what she is doing
1101              useFpump=action;
1102              break;
1103            case CUTSSTRATEGY:
1104              gomoryAction = action;
1105              probingAction = action;
1106              knapsackAction = action;
1107              cliqueAction = action;
1108              flowAction = action;
1109              mixedAction = action;
1110              twomirAction = action;
1111              parameters[whichParam(GOMORYCUTS,numberParameters,parameters)].setCurrentOption(action);
1112              parameters[whichParam(PROBINGCUTS,numberParameters,parameters)].setCurrentOption(action);
1113              parameters[whichParam(KNAPSACKCUTS,numberParameters,parameters)].setCurrentOption(action);
1114              if (!action) {
1115                redsplitAction = action;
1116                parameters[whichParam(REDSPLITCUTS,numberParameters,parameters)].setCurrentOption(action);
1117              }
1118              parameters[whichParam(CLIQUECUTS,numberParameters,parameters)].setCurrentOption(action);
1119              parameters[whichParam(FLOWCUTS,numberParameters,parameters)].setCurrentOption(action);
1120              parameters[whichParam(MIXEDCUTS,numberParameters,parameters)].setCurrentOption(action);
1121              parameters[whichParam(TWOMIRCUTS,numberParameters,parameters)].setCurrentOption(action);
1122              break;
1123            case HEURISTICSTRATEGY:
1124              useRounding = action;
1125              useGreedy = action;
1126              useCombine = action;
1127              //useLocalTree = action;
1128              useFpump=action;
1129              parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption(action);
1130              parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption(action);
1131              parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption(action);
1132              //parameters[whichParam(LOCALTREE,numberParameters,parameters)].setCurrentOption(action);
1133              parameters[whichParam(FPUMP,numberParameters,parameters)].setCurrentOption(action);
1134              break;
1135            case GREEDY:
1136              defaultSettings=false; // user knows what she is doing
1137              useGreedy = action;
1138              break;
1139            case COMBINE:
1140              defaultSettings=false; // user knows what she is doing
1141              useCombine = action;
1142              break;
1143            case LOCALTREE:
1144              defaultSettings=false; // user knows what she is doing
1145              useLocalTree = action;
1146              break;
1147            case COSTSTRATEGY:
1148              if (action!=1&&action!=0) {
1149                printf("Pseudo costs not implemented yet\n");
1150              } else {
1151                useCosts=action;
1152              }
1153              break;
1154            case PREPROCESS:
1155              preProcess = action;
1156              break;
1157            default:
1158              abort();
1159            }
1160          }
1161        } else {
1162          // action
1163          if (type==EXIT) {
1164#ifdef CBC_AMPL
1165            if(usingAmpl) {
1166              if (info.numberIntegers||info.numberBinary) {
1167                // integer
1168              } else {
1169                // linear
1170              }
1171              writeAmpl(&info);
1172              freeArrays2(&info);
1173              freeArgs(&info);
1174            }
1175#endif
1176            break; // stop all
1177          }
1178          switch (type) {
1179          case DUALSIMPLEX:
1180          case PRIMALSIMPLEX:
1181          case SOLVECONTINUOUS:
1182          case BARRIER:
1183            if (goodModel) {
1184              double objScale = 
1185                parameters[whichParam(OBJSCALE2,numberParameters,parameters)].doubleValue();
1186              if (objScale!=1.0) {
1187                int iColumn;
1188                int numberColumns=lpSolver->numberColumns();
1189                double * dualColumnSolution = 
1190                  lpSolver->dualColumnSolution();
1191                ClpObjective * obj = lpSolver->objectiveAsObject();
1192                assert(dynamic_cast<ClpLinearObjective *> (obj));
1193                double offset;
1194                double * objective = obj->gradient(NULL,NULL,offset,true);
1195                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1196                  dualColumnSolution[iColumn] *= objScale;
1197                  objective[iColumn] *= objScale;;
1198                }
1199                int iRow;
1200                int numberRows=lpSolver->numberRows();
1201                double * dualRowSolution = 
1202                  lpSolver->dualRowSolution();
1203                for (iRow=0;iRow<numberRows;iRow++) 
1204                  dualRowSolution[iRow] *= objScale;
1205                lpSolver->setObjectiveOffset(objScale*lpSolver->objectiveOffset());
1206              }
1207              ClpSolve::SolveType method;
1208              ClpSolve::PresolveType presolveType;
1209              ClpSimplex * model2 = lpSolver;
1210              if (dualize) {
1211                model2 = ((ClpSimplexOther *) model2)->dualOfModel();
1212                printf("Dual of model has %d rows and %d columns\n",
1213                       model2->numberRows(),model2->numberColumns());
1214                model2->setOptimizationDirection(1.0);
1215              }
1216              if (noPrinting)
1217                lpSolver->setLogLevel(0);
1218              ClpSolve solveOptions;
1219              solveOptions.setPresolveActions(presolveOptions);
1220              solveOptions.setSubstitution(substitution);
1221              if (preSolve!=5&&preSolve) {
1222                presolveType=ClpSolve::presolveNumber;
1223                if (preSolve<0) {
1224                  preSolve = - preSolve;
1225                  if (preSolve<=100) {
1226                    presolveType=ClpSolve::presolveNumber;
1227                    printf("Doing %d presolve passes - picking up non-costed slacks\n",
1228                           preSolve);
1229                    solveOptions.setDoSingletonColumn(true);
1230                  } else {
1231                    preSolve -=100;
1232                    presolveType=ClpSolve::presolveNumberCost;
1233                    printf("Doing %d presolve passes - picking up costed slacks\n",
1234                           preSolve);
1235                  }
1236                } 
1237              } else if (preSolve) {
1238                presolveType=ClpSolve::presolveOn;
1239              } else {
1240                presolveType=ClpSolve::presolveOff;
1241              }
1242              solveOptions.setPresolveType(presolveType,preSolve);
1243              if (type==DUALSIMPLEX||type==SOLVECONTINUOUS) {
1244                method=ClpSolve::useDual;
1245              } else if (type==PRIMALSIMPLEX) {
1246                method=ClpSolve::usePrimalorSprint;
1247              } else {
1248                method = ClpSolve::useBarrier;
1249                if (crossover==1) {
1250                  method=ClpSolve::useBarrierNoCross;
1251                } else if (crossover==2) {
1252                  ClpObjective * obj = lpSolver->objectiveAsObject();
1253                  if (obj->type()>1) {
1254                    method=ClpSolve::useBarrierNoCross;
1255                    presolveType=ClpSolve::presolveOff;
1256                    solveOptions.setPresolveType(presolveType,preSolve);
1257                  } 
1258                }
1259              }
1260              solveOptions.setSolveType(method);
1261              if(preSolveFile)
1262                presolveOptions |= 0x40000000;
1263              solveOptions.setSpecialOption(4,presolveOptions);
1264              solveOptions.setSpecialOption(5,printOptions);
1265              if (method==ClpSolve::useDual) {
1266                // dual
1267                if (doCrash)
1268                  solveOptions.setSpecialOption(0,1,doCrash); // crash
1269                else if (doIdiot)
1270                  solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
1271              } else if (method==ClpSolve::usePrimalorSprint) {
1272                // primal
1273                // if slp turn everything off
1274                if (slpValue>0) {
1275                  doCrash=false;
1276                  doSprint=0;
1277                  doIdiot=-1;
1278                  solveOptions.setSpecialOption(1,10,slpValue); // slp
1279                  method=ClpSolve::usePrimal;
1280                }
1281                if (doCrash) {
1282                  solveOptions.setSpecialOption(1,1,doCrash); // crash
1283                } else if (doSprint>0) {
1284                  // sprint overrides idiot
1285                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
1286                } else if (doIdiot>0) {
1287                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
1288                } else if (slpValue<=0) {
1289                  if (doIdiot==0) {
1290                    if (doSprint==0)
1291                      solveOptions.setSpecialOption(1,4); // all slack
1292                    else
1293                      solveOptions.setSpecialOption(1,9); // all slack or sprint
1294                  } else {
1295                    if (doSprint==0)
1296                      solveOptions.setSpecialOption(1,8); // all slack or idiot
1297                    else
1298                      solveOptions.setSpecialOption(1,7); // initiative
1299                  }
1300                }
1301                if (basisHasValues==-1)
1302                  solveOptions.setSpecialOption(1,11); // switch off values
1303              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
1304                int barrierOptions = choleskyType;
1305                if (scaleBarrier)
1306                  barrierOptions |= 8;
1307                if (doKKT)
1308                  barrierOptions |= 16;
1309                if (gamma)
1310                  barrierOptions |= 32*gamma;
1311                if (crossover==3) 
1312                  barrierOptions |= 256; // try presolve in crossover
1313                solveOptions.setSpecialOption(4,barrierOptions);
1314              }
1315              model2->initialSolve(solveOptions);
1316              basisHasValues=1;
1317              if (dualize) {
1318                ((ClpSimplexOther *) lpSolver)->restoreFromDual(model2);
1319                delete model2;
1320                model2=lpSolver;
1321              }
1322#ifdef CBC_AMPL
1323              if (usingAmpl) {
1324                double value = model2->getObjValue()*model2->getObjSense();
1325                char buf[300];
1326                int pos=0;
1327                int iStat = model2->status();
1328                if (iStat==0) {
1329                  pos += sprintf(buf+pos,"optimal," );
1330                } else if (iStat==1) {
1331                  // infeasible
1332                  pos += sprintf(buf+pos,"infeasible,");
1333                } else if (iStat==2) {
1334                  // unbounded
1335                  pos += sprintf(buf+pos,"unbounded,");
1336                } else if (iStat==3) {
1337                  pos += sprintf(buf+pos,"stopped on iterations or time,");
1338                } else if (iStat==4) {
1339                  iStat = 7;
1340                  pos += sprintf(buf+pos,"stopped on difficulties,");
1341                } else if (iStat==5) {
1342                  iStat = 3;
1343                  pos += sprintf(buf+pos,"stopped on ctrl-c,");
1344                } else {
1345                  pos += sprintf(buf+pos,"status unknown,");
1346                  iStat=6;
1347                }
1348                info.problemStatus=iStat;
1349                info.objValue = value;
1350                pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
1351                               value);
1352                sprintf(buf+pos,"\n%d iterations",
1353                        model2->getIterationCount());
1354                free(info.primalSolution);
1355                int numberColumns=model2->numberColumns();
1356                info.primalSolution = (double *) malloc(numberColumns*sizeof(double));
1357                CoinCopyN(model2->primalColumnSolution(),numberColumns,info.primalSolution);
1358                int numberRows = model2->numberRows();
1359                free(info.dualSolution);
1360                info.dualSolution = (double *) malloc(numberRows*sizeof(double));
1361                CoinCopyN(model2->dualRowSolution(),numberRows,info.dualSolution);
1362                CoinWarmStartBasis * basis = model2->getBasis();
1363                free(info.rowStatus);
1364                info.rowStatus = (int *) malloc(numberRows*sizeof(int));
1365                free(info.columnStatus);
1366                info.columnStatus = (int *) malloc(numberColumns*sizeof(int));
1367                // Put basis in
1368                int i;
1369                // free,basic,ub,lb are 0,1,2,3
1370                for (i=0;i<numberRows;i++) {
1371                  CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
1372                  info.rowStatus[i]=status;
1373                }
1374                for (i=0;i<numberColumns;i++) {
1375                  CoinWarmStartBasis::Status status = basis->getStructStatus(i);
1376                  info.columnStatus[i]=status;
1377                }
1378                // put buffer into info
1379                strcpy(info.buffer,buf);
1380                delete basis;
1381              }
1382#endif
1383            } else {
1384              std::cout<<"** Current model not valid"<<std::endl;
1385            }
1386            break;
1387          case STATISTICS:
1388            if (goodModel) {
1389              // If presolve on look at presolved
1390              bool deleteModel2=false;
1391              ClpSimplex * model2 = lpSolver;
1392              if (preSolve) {
1393                ClpPresolve pinfo;
1394                int presolveOptions2 = presolveOptions&~0x40000000;
1395                if ((presolveOptions2&0xffff)!=0)
1396                  pinfo.setPresolveActions(presolveOptions2);
1397                pinfo.setSubstitution(substitution);
1398                if ((printOptions&1)!=0)
1399                  pinfo.statistics();
1400                double presolveTolerance = 
1401                  parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1402                model2 = 
1403                  pinfo.presolvedModel(*lpSolver,presolveTolerance,
1404                                       true,preSolve);
1405                if (model2) {
1406                  printf("Statistics for presolved model\n");
1407                  deleteModel2=true;
1408                } else {
1409                  printf("Presolved model looks infeasible - will use unpresolved\n");
1410                  model2 = lpSolver;
1411                }
1412              } else {
1413                printf("Statistics for unpresolved model\n");
1414                model2 =  lpSolver;
1415              }
1416              statistics(lpSolver,model2);
1417              if (deleteModel2)
1418                delete model2;
1419            } else {
1420              std::cout<<"** Current model not valid"<<std::endl;
1421            }
1422            break;
1423          case TIGHTEN:
1424            if (goodModel) {
1425              int numberInfeasibilities = lpSolver->tightenPrimalBounds();
1426              if (numberInfeasibilities)
1427                std::cout<<"** Analysis indicates model infeasible"<<std::endl;
1428            } else {
1429              std::cout<<"** Current model not valid"<<std::endl;
1430            }
1431            break;
1432          case PLUSMINUS:
1433            if (goodModel) {
1434              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
1435              ClpPackedMatrix* clpMatrix =
1436                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1437              if (clpMatrix) {
1438                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
1439                if (newMatrix->getIndices()) {
1440                  lpSolver->replaceMatrix(newMatrix);
1441                  delete saveMatrix;
1442                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
1443                } else {
1444                  std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
1445                }
1446              } else {
1447                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
1448              }
1449            } else {
1450              std::cout<<"** Current model not valid"<<std::endl;
1451            }
1452            break;
1453          case NETWORK:
1454            if (goodModel) {
1455              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
1456              ClpPackedMatrix* clpMatrix =
1457                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1458              if (clpMatrix) {
1459                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
1460                if (newMatrix->getIndices()) {
1461                  lpSolver->replaceMatrix(newMatrix);
1462                  delete saveMatrix;
1463                  std::cout<<"Matrix converted to network matrix"<<std::endl;
1464                } else {
1465                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
1466                }
1467              } else {
1468                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
1469              }
1470            } else {
1471              std::cout<<"** Current model not valid"<<std::endl;
1472            }
1473            break;
1474          case MIPLIB:
1475            // User can set options - main differenec is lack of model and CglPreProcess
1476            goodModel=true;
1477/*
1478  Run branch-and-cut. First set a few options -- node comparison, scaling. If
1479  the solver is Clp, consider running some presolve code (not yet converted
1480  this to generic OSI) with branch-and-cut. If presolve is disabled, or the
1481  solver is not Clp, simply run branch-and-cut. Print elapsed time at the end.
1482*/
1483          case BAB: // branchAndBound
1484          case STRENGTHEN:
1485            if (goodModel) {
1486              bool miplib = type==MIPLIB;
1487              int logLevel = parameters[slog].intValue();
1488              // Reduce printout
1489              if (logLevel<=1)
1490                model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
1491              // Don't switch off all output
1492              {
1493                OsiSolverInterface * solver = model.solver();
1494                OsiClpSolverInterface * si =
1495                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
1496                assert (si != NULL);
1497                si->setSpecialOptions(0x40000000);
1498              }
1499              if (!miplib)
1500                model.initialSolve();
1501              // If user made settings then use them
1502              if (!defaultSettings) {
1503                OsiSolverInterface * solver = model.solver();
1504                if (!doScaling)
1505                  solver->setHintParam(OsiDoScale,false,OsiHintTry);
1506                OsiClpSolverInterface * si =
1507                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
1508                assert (si != NULL);
1509                // get clp itself
1510                ClpSimplex * modelC = si->getModelPtr();
1511                //if (modelC->tightenPrimalBounds()!=0) {
1512                //std::cout<<"Problem is infeasible!"<<std::endl;
1513                //break;
1514                //}
1515                // bounds based on continuous
1516                if (tightenFactor) {
1517                  if (modelC->tightenPrimalBounds(tightenFactor)!=0) {
1518                    std::cout<<"Problem is infeasible!"<<std::endl;
1519                    break;
1520                  }
1521                }
1522                if (djFix<1.0e20) {
1523                  // do some fixing
1524                  int numberColumns = modelC->numberColumns();
1525                  int i;
1526                  const char * type = modelC->integerInformation();
1527                  double * lower = modelC->columnLower();
1528                  double * upper = modelC->columnUpper();
1529                  double * solution = modelC->primalColumnSolution();
1530                  double * dj = modelC->dualColumnSolution();
1531                  int numberFixed=0;
1532                  for (i=0;i<numberColumns;i++) {
1533                    if (type[i]) {
1534                      double value = solution[i];
1535                      if (value<lower[i]+1.0e-5&&dj[i]>djFix) {
1536                        solution[i]=lower[i];
1537                        upper[i]=lower[i];
1538                        numberFixed++;
1539                      } else if (value>upper[i]-1.0e-5&&dj[i]<-djFix) {
1540                        solution[i]=upper[i];
1541                        lower[i]=upper[i];
1542                        numberFixed++;
1543                      }
1544                    }
1545                  }
1546                  printf("%d columns fixed\n",numberFixed);
1547                }
1548              }
1549              // See if we want preprocessing
1550              OsiSolverInterface * saveSolver=NULL;
1551              CglPreProcess process;
1552              delete babModel;
1553              babModel = new CbcModel(model);
1554              OsiSolverInterface * solver3 = clpSolver->clone();
1555              babModel->assignSolver(solver3);
1556              OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
1557              int numberChanged=0;
1558              if (clpSolver2->messageHandler()->logLevel())
1559                clpSolver2->messageHandler()->setLogLevel(1);
1560              if (logLevel>-1)
1561                clpSolver2->messageHandler()->setLogLevel(logLevel);
1562              lpSolver = clpSolver2->getModelPtr();
1563              if (lpSolver->factorizationFrequency()==200&&!miplib) {
1564                // User did not touch preset
1565                int numberRows = lpSolver->numberRows();
1566                const int cutoff1=10000;
1567                const int cutoff2=100000;
1568                const int base=75;
1569                const int freq0 = 50;
1570                const int freq1=200;
1571                const int freq2=400;
1572                const int maximum=1000;
1573                int frequency;
1574                if (numberRows<cutoff1)
1575                  frequency=base+numberRows/freq0;
1576                else if (numberRows<cutoff2)
1577                  frequency=base+cutoff1/freq0 + (numberRows-cutoff1)/freq1;
1578                else
1579                  frequency=base+cutoff1/freq0 + (cutoff2-cutoff1)/freq1 + (numberRows-cutoff2)/freq2;
1580                lpSolver->setFactorizationFrequency(CoinMin(maximum,frequency));
1581              }
1582              time2 = CoinCpuTime();
1583              totalTime += time2-time1;
1584              time1 = time2;
1585              double timeLeft = babModel->getMaximumSeconds();
1586              if (preProcess&&type==BAB) {
1587                saveSolver=babModel->solver()->clone();
1588                /* Do not try and produce equality cliques and
1589                   do up to 10 passes */
1590                OsiSolverInterface * solver2;
1591                {
1592                  // Tell solver we are in Branch and Cut
1593                  saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;
1594                  // Default set of cut generators
1595                  CglProbing generator1;
1596                  generator1.setUsingObjective(true);
1597                  generator1.setMaxPass(3);
1598                  generator1.setMaxProbeRoot(saveSolver->getNumCols());
1599                  generator1.setMaxElements(100);
1600                  generator1.setMaxLookRoot(50);
1601                  generator1.setRowCuts(3);
1602                  // Add in generators
1603                  process.addCutGenerator(&generator1);
1604                  int translate[]={9999,0,0,-1,2,3};
1605                  process.messageHandler()->setLogLevel(babModel->logLevel());
1606                  solver2 = process.preProcessNonDefault(*saveSolver,translate[preProcess],10);
1607                  // Tell solver we are not in Branch and Cut
1608                  saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
1609                  if (solver2)
1610                    solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
1611                }
1612#ifdef CBC_AMPL
1613                if (!solver2&&usingAmpl) {
1614                  // infeasible
1615                  info.problemStatus=1;
1616                  info.objValue = 1.0e100;
1617                  sprintf(info.buffer,"infeasible by pre-processing");
1618                  info.primalSolution=NULL;
1619                  info.dualSolution=NULL;
1620                  break;
1621                }
1622#endif
1623                if (!noPrinting) {
1624                  if (!solver2) {
1625                    printf("Pre-processing says infeasible\n");
1626                    break;
1627                  } else {
1628                    printf("processed model has %d rows, %d columns and %d elements\n",
1629                           solver2->getNumRows(),solver2->getNumCols(),solver2->getNumElements());
1630                  }
1631                }
1632                //solver2->resolve();
1633                if (preProcess==2) {
1634                  OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (solver2);
1635                  ClpSimplex * lpSolver = clpSolver2->getModelPtr();
1636                  lpSolver->writeMps("presolved.mps",0,1,lpSolver->optimizationDirection());
1637                  printf("Preprocessed model (minimization) on presolved.mps\n");
1638                }
1639                // we have to keep solver2 so pass clone
1640                solver2 = solver2->clone();
1641                babModel->assignSolver(solver2);
1642                babModel->initialSolve();
1643                babModel->setMaximumSeconds(timeLeft-(CoinCpuTime()-time1));
1644              }
1645              // now tighten bounds
1646              if (!miplib) {
1647                OsiClpSolverInterface * si =
1648                  dynamic_cast<OsiClpSolverInterface *>(babModel->solver()) ;
1649                assert (si != NULL);
1650                // get clp itself
1651                ClpSimplex * modelC = si->getModelPtr();
1652                if (noPrinting)
1653                  modelC->setLogLevel(0);
1654                if (modelC->tightenPrimalBounds()!=0) {
1655                  std::cout<<"Problem is infeasible!"<<std::endl;
1656                  break;
1657                }
1658                modelC->dual();
1659              }
1660              if (debugValues) {
1661                // for debug
1662                std::string problemName ;
1663                babModel->solver()->getStrParam(OsiProbName,problemName) ;
1664                //babModel->solver()->activateRowCutDebugger(problemName.c_str()) ;
1665                twomirGen.probname_=strdup(problemName.c_str());
1666                // checking seems odd
1667                //redsplitGen.set_given_optsol(babModel->solver()->getRowCutDebuggerAlways()->optimalSolution(),
1668                //                         babModel->getNumCols());
1669              }
1670              if (useCosts) {
1671                int numberColumns = babModel->getNumCols();
1672                int * sort = new int[numberColumns];
1673                double * dsort = new double[numberColumns];
1674                int * priority = new int [numberColumns];
1675                const double * objective = babModel->getObjCoefficients();
1676                int iColumn;
1677                int n=0;
1678                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1679                  if (babModel->isInteger(iColumn)) {
1680                    sort[n]=n;
1681                    dsort[n++]=-objective[iColumn];
1682                  }
1683                }
1684                CoinSort_2(dsort,dsort+n,sort);
1685                int level=0;
1686                double last = -1.0e100;
1687                for (int i=0;i<n;i++) {
1688                  int iPut=sort[i];
1689                  if (dsort[i]!=last) {
1690                    level++;
1691                    last=dsort[i];
1692                  }
1693                  priority[iPut]=level;
1694                }
1695                babModel->passInPriorities( priority,false);
1696                delete [] priority;
1697                delete [] sort;
1698                delete [] dsort;
1699              }
1700              // FPump done first as it only works if no solution
1701              CbcHeuristicFPump heuristic4(*babModel);
1702              if (useFpump) {
1703                heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
1704                babModel->addHeuristic(&heuristic4);
1705              }
1706              if (!miplib) {
1707                CbcRounding heuristic1(*babModel);
1708                if (useRounding)
1709                  babModel->addHeuristic(&heuristic1) ;
1710                CbcHeuristicLocal heuristic2(*babModel);
1711                heuristic2.setSearchType(1);
1712                if (useCombine)
1713                  babModel->addHeuristic(&heuristic2);
1714                CbcHeuristicGreedyCover heuristic3(*babModel);
1715                CbcHeuristicGreedyEquality heuristic3a(*babModel);
1716                if (useGreedy) {
1717                  babModel->addHeuristic(&heuristic3);
1718                  babModel->addHeuristic(&heuristic3a);
1719                }
1720                if (useLocalTree) {
1721                  CbcTreeLocal localTree(babModel,NULL,10,0,0,10000,2000);
1722                  babModel->passInTreeHandler(localTree);
1723                }
1724              }
1725              // add cut generators if wanted
1726              int switches[20];
1727              int numberGenerators=0;
1728              if (probingAction==1) {
1729                babModel->addCutGenerator(&probingGen,-1,"Probing");
1730                switches[numberGenerators++]=0;
1731              } else if (probingAction>=2) {
1732                babModel->addCutGenerator(&probingGen,-101+probingAction,"Probing");
1733                switches[numberGenerators++]=0;
1734              }
1735              if (gomoryAction==1) {
1736                babModel->addCutGenerator(&gomoryGen,-1,"Gomory");
1737                switches[numberGenerators++]=1;
1738              } else if (gomoryAction>=2) {
1739                babModel->addCutGenerator(&gomoryGen,-101+gomoryAction,"Gomory");
1740                switches[numberGenerators++]=-1;
1741              }
1742              if (knapsackAction==1) {
1743                babModel->addCutGenerator(&knapsackGen,-1,"Knapsack");
1744                switches[numberGenerators++]=0;
1745              } else if (knapsackAction>=2) {
1746                babModel->addCutGenerator(&knapsackGen,-101+knapsackAction,"Knapsack");
1747                switches[numberGenerators++]=0;
1748              }
1749              if (redsplitAction==1) {
1750                babModel->addCutGenerator(&redsplitGen,-1,"Reduce-and-split");
1751                switches[numberGenerators++]=1;
1752              } else if (redsplitAction>=2) {
1753                babModel->addCutGenerator(&redsplitGen,-101+redsplitAction,"Reduce-and-split");
1754                switches[numberGenerators++]=1;
1755              }
1756              if (cliqueAction==1) {
1757                babModel->addCutGenerator(&cliqueGen,-1,"Clique");
1758                switches[numberGenerators++]=1;
1759              } else if (cliqueAction>=2) {
1760                babModel->addCutGenerator(&cliqueGen,-101+cliqueAction,"Clique");
1761                switches[numberGenerators++]=-1;
1762              }
1763              if (mixedAction==1) {
1764                babModel->addCutGenerator(&mixedGen,-1,"MixedIntegerRounding2");
1765                switches[numberGenerators++]=1;
1766              } else if (mixedAction>=2) {
1767                babModel->addCutGenerator(&mixedGen,-101+mixedAction,"MixedIntegerRounding2");
1768                switches[numberGenerators++]=-1;
1769              }
1770              if (flowAction==1) {
1771                babModel->addCutGenerator(&flowGen,-1,"FlowCover");
1772                switches[numberGenerators++]=1;
1773              } else if (flowAction>=2) {
1774                babModel->addCutGenerator(&flowGen,-101+flowAction,"FlowCover");
1775                switches[numberGenerators++]=1;
1776              }
1777              if (twomirAction==1) {
1778                babModel->addCutGenerator(&twomirGen,-1,"TwoMirCuts");
1779                switches[numberGenerators++]=1;
1780              } else if (twomirAction>=2) {
1781                babModel->addCutGenerator(&twomirGen,-101+twomirAction,"TwoMirCuts");
1782                switches[numberGenerators++]=1;
1783              }
1784              // Say we want timings
1785              numberGenerators = babModel->numberCutGenerators();
1786              int iGenerator;
1787              int cutDepth=
1788                parameters[whichParam(CUTDEPTH,numberParameters,parameters)].intValue();
1789              for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
1790                CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
1791                int howOften = generator->howOften();
1792                if (howOften==-98||howOften==-99) 
1793                  generator->setSwitchOffIfLessThan(switches[iGenerator]);
1794                generator->setTiming(true);
1795                if (cutDepth>=0)
1796                  generator->setWhatDepth(cutDepth) ;
1797              }
1798              // Could tune more
1799              if (!miplib) {
1800                babModel->setMinimumDrop(min(5.0e-2,
1801                                             fabs(babModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
1802                if (cutPass==-1234567) {
1803                  if (babModel->getNumCols()<500)
1804                    babModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
1805                  else if (babModel->getNumCols()<5000)
1806                    babModel->setMaximumCutPassesAtRoot(100); // use minimum drop
1807                  else
1808                    babModel->setMaximumCutPassesAtRoot(20);
1809                } else {
1810                  babModel->setMaximumCutPassesAtRoot(cutPass);
1811                }
1812                babModel->setMaximumCutPasses(1);
1813              }
1814              // Do more strong branching if small
1815              //if (babModel->getNumCols()<5000)
1816              //babModel->setNumberStrong(20);
1817              // Switch off strong branching if wanted
1818              //if (babModel->getNumCols()>10*babModel->getNumRows())
1819              //babModel->setNumberStrong(0);
1820              if (!noPrinting) {
1821                babModel->messageHandler()->setLogLevel(parameters[log].intValue());
1822                if (babModel->getNumCols()>2000||babModel->getNumRows()>1500||
1823                    babModel->messageHandler()->logLevel()>1)
1824                  babModel->setPrintFrequency(100);
1825              }
1826             
1827              babModel->solver()->setIntParam(OsiMaxNumIterationHotStart,
1828                    parameters[whichParam(MAXHOTITS,numberParameters,parameters)].intValue());
1829              OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
1830              // go faster stripes
1831              if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
1832                osiclp->setupForRepeatedUse(2,parameters[slog].intValue());
1833              } else {
1834                osiclp->setupForRepeatedUse(0,parameters[slog].intValue());
1835              }
1836              double increment=babModel->getCutoffIncrement();;
1837              int * changed = NULL;
1838              if (!miplib)
1839                changed=analyze( osiclp,numberChanged,increment,false);
1840              if (debugValues) {
1841                if (numberDebugValues==babModel->getNumCols()) {
1842                  // for debug
1843                  babModel->solver()->activateRowCutDebugger(debugValues) ;
1844                } else {
1845                  printf("debug file has incorrect number of columns\n");
1846                }
1847              }
1848              babModel->setCutoffIncrement(CoinMax(babModel->getCutoffIncrement(),increment));
1849              // Turn this off if you get problems
1850              // Used to be automatically set
1851              int mipOptions = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue();
1852              if (mipOptions!=(128|64|1))
1853                printf("mip options %d\n",mipOptions);
1854              osiclp->setSpecialOptions(mipOptions);
1855              if (gapRatio < 1.0e100) {
1856                double value = babModel->solver()->getObjValue() ;
1857                double value2 = gapRatio*(1.0e-5+fabs(value)) ;
1858                babModel->setAllowableGap(value2) ;
1859                std::cout << "Continuous " << value
1860                          << ", so allowable gap set to "
1861                          << value2 << std::endl ;
1862              }
1863              // probably faster to use a basis to get integer solutions
1864              babModel->setSpecialOptions(2);
1865              currentBranchModel = babModel;
1866              OsiSolverInterface * strengthenedModel=NULL;
1867              if (type==BAB||type==MIPLIB) {
1868                int moreMipOptions = parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].intValue();
1869                if (moreMipOptions>=0) {
1870                  printf("more mip options %d\n",moreMipOptions);
1871                  babModel->setSearchStrategy(moreMipOptions);
1872                }
1873              }
1874              if (type==BAB) {
1875#ifdef CBC_AMPL
1876                if (usingAmpl) {
1877                  priorities=info.priorities;
1878                  branchDirection=info.branchDirection;
1879                  pseudoDown=info.pseudoDown;
1880                  pseudoUp=info.pseudoUp;
1881                }
1882#endif               
1883                const int * originalColumns = NULL;
1884                if (preProcess)
1885                  originalColumns = process.originalColumns();
1886                if (preProcess&&process.numberSOS()) {
1887                  int numberSOS = process.numberSOS();
1888                  int numberIntegers = babModel->numberIntegers();
1889                  /* model may not have created objects
1890                     If none then create
1891                  */
1892                  if (!numberIntegers||!babModel->numberObjects()) {
1893                    int type = (pseudoUp) ? 1 : 0;
1894                    babModel->findIntegers(true,type);
1895                    numberIntegers = babModel->numberIntegers();
1896                  }
1897                  CbcObject ** oldObjects = babModel->objects();
1898                  // Do sets and priorities
1899                  CbcObject ** objects = new CbcObject * [numberSOS];
1900                  // set old objects to have low priority
1901                  int numberOldObjects = babModel->numberObjects();
1902                  int numberColumns = babModel->getNumCols();
1903                  for (int iObj = 0;iObj<numberOldObjects;iObj++) {
1904                    oldObjects[iObj]->setPriority(numberColumns+1);
1905                    int iColumn = oldObjects[iObj]->columnNumber();
1906                    assert (iColumn>=0);
1907                    if (originalColumns)
1908                      iColumn = originalColumns[iColumn];
1909                    if (branchDirection)
1910                      oldObjects[iObj]->setPreferredWay(branchDirection[iColumn]);
1911                    if (pseudoUp) {
1912                      CbcSimpleIntegerPseudoCost * obj1a =
1913                        dynamic_cast <CbcSimpleIntegerPseudoCost *>(oldObjects[iObj]) ;
1914                      assert (obj1a);
1915                      obj1a->setDownPseudoCost(pseudoDown[iColumn]);
1916                      obj1a->setUpPseudoCost(pseudoUp[iColumn]);
1917                    }
1918                  }
1919                  const int * starts = process.startSOS();
1920                  const int * which = process.whichSOS();
1921                  const int * type = process.typeSOS();
1922                  const double * weight = process.weightSOS();
1923                  int iSOS;
1924                  for (iSOS =0;iSOS<numberSOS;iSOS++) {
1925                    int iStart = starts[iSOS];
1926                    int n=starts[iSOS+1]-iStart;
1927                    objects[iSOS] = new CbcSOS(babModel,n,which+iStart,weight+iStart,
1928                                               iSOS,type[iSOS]);
1929                    // branch on long sets first
1930                    objects[iSOS]->setPriority(numberColumns-n);
1931                  }
1932                  babModel->addObjects(numberSOS,objects);
1933                  for (iSOS=0;iSOS<numberSOS;iSOS++)
1934                    delete objects[iSOS];
1935                  delete [] objects;
1936                } else if (priorities||branchDirection||pseudoDown||pseudoUp) {
1937                  // do anyway for priorities etc
1938                  int numberIntegers = babModel->numberIntegers();
1939                  /* model may not have created objects
1940                     If none then create
1941                  */
1942                  if (!numberIntegers||!babModel->numberObjects()) {
1943                    int type = (pseudoUp) ? 1 : 0;
1944                    babModel->findIntegers(true,type);
1945                  }
1946                  CbcObject ** objects = babModel->objects();
1947                  int numberObjects = babModel->numberObjects();
1948                  for (int iObj = 0;iObj<numberObjects;iObj++) {
1949                    int iColumn = objects[iObj]->columnNumber();
1950                    assert (iColumn>=0);
1951                    if (originalColumns)
1952                      iColumn = originalColumns[iColumn];
1953                    if (branchDirection)
1954                      objects[iObj]->setPreferredWay(branchDirection[iColumn]);
1955                    if (priorities) {
1956                      int iPriority = priorities[iColumn];
1957                      if (iPriority>0)
1958                        objects[iObj]->setPriority(iPriority);
1959                    }
1960                    if (pseudoUp&&pseudoUp[iColumn]) {
1961                      CbcSimpleIntegerPseudoCost * obj1a =
1962                        dynamic_cast <CbcSimpleIntegerPseudoCost *>(objects[iObj]) ;
1963                      assert (obj1a);
1964                      obj1a->setDownPseudoCost(pseudoDown[iColumn]);
1965                      obj1a->setUpPseudoCost(pseudoUp[iColumn]);
1966                    }
1967                  }
1968                }
1969#ifdef CBC_AMPL
1970                if (!usingAmpl) {
1971#endif
1972                  free(priorities);
1973                  priorities=NULL;
1974                  free(branchDirection);
1975                  branchDirection=NULL;
1976                  free(pseudoDown);
1977                  pseudoDown=NULL;
1978                  free(pseudoUp);
1979                  pseudoUp=NULL;
1980#ifdef CBC_AMPL
1981                }
1982#endif               
1983                int statistics = (printOptions>0) ? printOptions: 0;
1984                babModel->branchAndBound(statistics);
1985              } else if (type==MIPLIB) {
1986                CbcStrategyDefault strategy(true,5,5);
1987                // Set up pre-processing to find sos if wanted
1988                if (preProcess)
1989                  strategy.setupPreProcessing(2);
1990                babModel->setStrategy(strategy);
1991                CbcClpUnitTest(*babModel);
1992                goodModel=false;
1993                break;
1994              } else {
1995                strengthenedModel = babModel->strengthenedModel();
1996              }
1997              currentBranchModel = NULL;
1998              osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
1999              if (debugFile=="createAfterPre"&&babModel->bestSolution()) {
2000                lpSolver = osiclp->getModelPtr();
2001                //move best solution (should be there -- but ..)
2002                int n = lpSolver->getNumCols();
2003                memcpy(lpSolver->primalColumnSolution(),babModel->bestSolution(),n*sizeof(double));
2004                saveSolution(osiclp->getModelPtr(),"debug.file");
2005              }
2006              if (!noPrinting) {
2007                // Print more statistics
2008                std::cout<<"Cuts at root node changed objective from "<<babModel->getContinuousObjective()
2009                         <<" to "<<babModel->rootObjectiveAfterCuts()<<std::endl;
2010               
2011                for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
2012                  CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
2013                  std::cout<<generator->cutGeneratorName()<<" was tried "
2014                           <<generator->numberTimesEntered()<<" times and created "
2015                           <<generator->numberCutsInTotal()<<" cuts of which "
2016                           <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
2017                  if (generator->timing())
2018                    std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
2019                  else
2020                    std::cout<<std::endl;
2021                }
2022              }
2023              time2 = CoinCpuTime();
2024              totalTime += time2-time1;
2025              // For best solution
2026              double * bestSolution = NULL;
2027              if (babModel->getMinimizationObjValue()<1.0e50&&type==BAB) {
2028                // post process
2029                if (preProcess) {
2030                  int n = saveSolver->getNumCols();
2031                  bestSolution = new double [n];
2032                  process.postProcess(*babModel->solver());
2033                  // Solution now back in saveSolver
2034                  babModel->assignSolver(saveSolver);
2035                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
2036                } else {
2037                  int n = babModel->solver()->getNumCols();
2038                  bestSolution = new double [n];
2039                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
2040                }
2041              }
2042              if (type==STRENGTHEN&&strengthenedModel)
2043                clpSolver = dynamic_cast< OsiClpSolverInterface*> (strengthenedModel);
2044              lpSolver = clpSolver->getModelPtr();
2045              if (numberChanged) {
2046                for (int i=0;i<numberChanged;i++) {
2047                  int iColumn=changed[i];
2048                  clpSolver->setContinuous(iColumn);
2049                }
2050                delete [] changed;
2051              }
2052              if (type==BAB) {
2053                //move best solution (should be there -- but ..)
2054                int n = lpSolver->getNumCols();
2055                if (bestSolution)
2056                  memcpy(lpSolver->primalColumnSolution(),bestSolution,n*sizeof(double));
2057                if (debugFile=="create"&&bestSolution) {
2058                  saveSolution(lpSolver,"debug.file");
2059                }
2060                delete [] bestSolution;
2061                std::string statusName[]={"Finished","Stopped on ","Difficulties",
2062                                          "","","User ctrl-c"};
2063                std::string minor[]={"","","gap","nodes","time","","solutions","user ctrl-c"};
2064                int iStat = babModel->status();
2065                int iStat2 = babModel->secondaryStatus();
2066                if (!noPrinting)
2067                  std::cout<<"Result - "<<statusName[iStat]<<minor[iStat2]
2068                           <<" objective "<<babModel->getObjValue()<<
2069                    " after "<<babModel->getNodeCount()<<" nodes and "
2070                           <<babModel->getIterationCount()<<
2071                    " iterations - took "<<time2-time1<<" seconds"<<std::endl;
2072#ifdef CBC_AMPL
2073                if (usingAmpl) {
2074                  double value = babModel->getObjValue()*lpSolver->getObjSense();
2075                  char buf[300];
2076                  int pos=0;
2077                  if (iStat==0) {
2078                    if (babModel->getObjValue()<1.0e40) {
2079                      pos += sprintf(buf+pos,"optimal," );
2080                    } else {
2081                      // infeasible
2082                      iStat=1;
2083                      pos += sprintf(buf+pos,"infeasible,");
2084                    }
2085                  } else if (iStat==1) {
2086                    if (iStat2!=6)
2087                      iStat=3;
2088                    else
2089                      iStat=4;
2090                    pos += sprintf(buf+pos,"stopped on %s,",minor[iStat2].c_str());
2091                  } else if (iStat==2) {
2092                    iStat = 7;
2093                    pos += sprintf(buf+pos,"stopped on difficulties,");
2094                  } else if (iStat==5) {
2095                    iStat = 3;
2096                    pos += sprintf(buf+pos,"stopped on ctrl-c,");
2097                  } else {
2098                    pos += sprintf(buf+pos,"status unknown,");
2099                    iStat=6;
2100                  }
2101                  info.problemStatus=iStat;
2102                  info.objValue = value;
2103                  if (babModel->getObjValue()<1.0e40) 
2104                    pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
2105                                   value);
2106                  sprintf(buf+pos,"\n%d nodes, %d iterations",
2107                          babModel->getNodeCount(),
2108                          babModel->getIterationCount());
2109                  if (bestSolution) {
2110                    free(info.primalSolution);
2111                    info.primalSolution = (double *) malloc(n*sizeof(double));
2112                    CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution);
2113                    int numberRows = lpSolver->numberRows();
2114                    free(info.dualSolution);
2115                    info.dualSolution = (double *) malloc(numberRows*sizeof(double));
2116                    CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution);
2117                  } else {
2118                    info.primalSolution=NULL;
2119                    info.dualSolution=NULL;
2120                  }
2121                  // put buffer into info
2122                  strcpy(info.buffer,buf);
2123                }
2124#endif
2125              } else {
2126                std::cout<<"Model strengthend - now has "<<clpSolver->getNumRows()
2127                         <<" rows"<<std::endl;
2128              }
2129              time1 = time2;
2130              delete babModel;
2131              babModel=NULL;
2132            } else {
2133              std::cout << "** Current model not valid" << std::endl ; 
2134            }
2135            break ;
2136          case IMPORT:
2137            {
2138#ifdef CBC_AMPL
2139              if (!usingAmpl) {
2140#endif
2141                free(priorities);
2142                priorities=NULL;
2143                free(branchDirection);
2144                branchDirection=NULL;
2145                free(pseudoDown);
2146                pseudoDown=NULL;
2147                free(pseudoUp);
2148                pseudoUp=NULL;
2149#ifdef CBC_AMPL
2150              }
2151#endif               
2152              delete babModel;
2153              babModel=NULL;
2154              // get next field
2155              field = CoinReadGetString(argc,argv);
2156              if (field=="$") {
2157                field = parameters[iParam].stringValue();
2158              } else if (field=="EOL") {
2159                parameters[iParam].printString();
2160                break;
2161              } else {
2162                parameters[iParam].setStringValue(field);
2163              }
2164              std::string fileName;
2165              bool canOpen=false;
2166              if (field=="-") {
2167                // stdin
2168                canOpen=true;
2169                fileName = "-";
2170              } else {
2171                bool absolutePath;
2172                if (dirsep=='/') {
2173                  // non Windows (or cygwin)
2174                  absolutePath=(field[0]=='/');
2175                } else {
2176                  //Windows (non cycgwin)
2177                  absolutePath=(field[0]=='\\');
2178                  // but allow for :
2179                  if (strchr(field.c_str(),':'))
2180                    absolutePath=true;
2181                }
2182                if (absolutePath) {
2183                  fileName = field;
2184                } else if (field[0]=='~') {
2185                  char * environVar = getenv("HOME");
2186                  if (environVar) {
2187                    std::string home(environVar);
2188                    field=field.erase(0,1);
2189                    fileName = home+field;
2190                  } else {
2191                    fileName=field;
2192                  }
2193                } else {
2194                  fileName = directory+field;
2195                }
2196                FILE *fp=fopen(fileName.c_str(),"r");
2197                if (fp) {
2198                  // can open - lets go for it
2199                  fclose(fp);
2200                  canOpen=true;
2201                } else {
2202                  std::cout<<"Unable to open file "<<fileName<<std::endl;
2203                }
2204              }
2205              if (canOpen) {
2206                int status =lpSolver->readMps(fileName.c_str(),
2207                                                   keepImportNames!=0,
2208                                                   allowImportErrors!=0);
2209                if (!status||(status>0&&allowImportErrors)) {
2210                  if (keepImportNames) {
2211                    lengthName = lpSolver->lengthNames();
2212                    rowNames = *(lpSolver->rowNames());
2213                    columnNames = *(lpSolver->columnNames());
2214                  } else {
2215                    lengthName=0;
2216                  }
2217                  goodModel=true;
2218                  //Set integers in clpsolver
2219                  const char * info = lpSolver->integerInformation();
2220                  if (info) {
2221                    int numberColumns = lpSolver->numberColumns();
2222                    int i;
2223                    for (i=0;i<numberColumns;i++) {
2224                      if (info[i]) 
2225                        clpSolver->setInteger(i);
2226                    }
2227                  }
2228                  // sets to all slack (not necessary?)
2229                  lpSolver->createStatus();
2230                  time2 = CoinCpuTime();
2231                  totalTime += time2-time1;
2232                  time1=time2;
2233                  // Go to canned file if just input file
2234                  if (CbcOrClpRead_mode==2&&argc==2) {
2235                    // only if ends .mps
2236                    std::string::size_type loc = fileName.find(".mps") ;
2237                    if (loc != std::string::npos &&
2238                        fileName.length() == loc+3)
2239                    { fileName.replace(loc+1,3,"par") ;
2240                      FILE *fp=fopen(fileName.c_str(),"r");
2241                      if (fp) {
2242                        CbcOrClpReadCommand=fp; // Read from that file
2243                        CbcOrClpRead_mode=-1;
2244                      }
2245                    }
2246                  }
2247                } else {
2248                  // errors
2249                  std::cout<<"There were "<<status<<
2250                    " errors on input"<<std::endl;
2251                }
2252              }
2253            }
2254            break;
2255          case EXPORT:
2256            if (goodModel) {
2257              // get next field
2258              field = CoinReadGetString(argc,argv);
2259              if (field=="$") {
2260                field = parameters[iParam].stringValue();
2261              } else if (field=="EOL") {
2262                parameters[iParam].printString();
2263                break;
2264              } else {
2265                parameters[iParam].setStringValue(field);
2266              }
2267              std::string fileName;
2268              bool canOpen=false;
2269              if (field[0]=='/'||field[0]=='\\') {
2270                fileName = field;
2271              } else if (field[0]=='~') {
2272                char * environVar = getenv("HOME");
2273                if (environVar) {
2274                  std::string home(environVar);
2275                  field=field.erase(0,1);
2276                  fileName = home+field;
2277                } else {
2278                  fileName=field;
2279                }
2280              } else {
2281                fileName = directory+field;
2282              }
2283              FILE *fp=fopen(fileName.c_str(),"w");
2284              if (fp) {
2285                // can open - lets go for it
2286                fclose(fp);
2287                canOpen=true;
2288              } else {
2289                std::cout<<"Unable to open file "<<fileName<<std::endl;
2290              }
2291              if (canOpen) {
2292                // If presolve on then save presolved
2293                bool deleteModel2=false;
2294                ClpSimplex * model2 = lpSolver;
2295                if (preSolve) {
2296                  ClpPresolve pinfo;
2297                  int presolveOptions2 = presolveOptions&~0x40000000;
2298                  if ((presolveOptions2&0xffff)!=0)
2299                    pinfo.setPresolveActions(presolveOptions2);
2300                  if ((printOptions&1)!=0)
2301                    pinfo.statistics();
2302                  double presolveTolerance = 
2303                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
2304                  model2 = 
2305                    pinfo.presolvedModel(*lpSolver,presolveTolerance,
2306                                         true,preSolve);
2307                  if (model2) {
2308                    printf("Saving presolved model on %s\n",
2309                           fileName.c_str());
2310                    deleteModel2=true;
2311                  } else {
2312                    printf("Presolved model looks infeasible - saving original on %s\n",
2313                           fileName.c_str());
2314                    deleteModel2=false;
2315                    model2 = lpSolver;
2316
2317                  }
2318                } else {
2319                  printf("Saving model on %s\n",
2320                           fileName.c_str());
2321                }
2322#if 0
2323                // Convert names
2324                int iRow;
2325                int numberRows=model2->numberRows();
2326                int iColumn;
2327                int numberColumns=model2->numberColumns();
2328
2329                char ** rowNames = NULL;
2330                char ** columnNames = NULL;
2331                if (model2->lengthNames()) {
2332                  rowNames = new char * [numberRows];
2333                  for (iRow=0;iRow<numberRows;iRow++) {
2334                    rowNames[iRow] =
2335                      strdup(model2->rowName(iRow).c_str());
2336#ifdef STRIPBLANKS
2337                    char * xx = rowNames[iRow];
2338                    int i;
2339                    int length = strlen(xx);
2340                    int n=0;
2341                    for (i=0;i<length;i++) {
2342                      if (xx[i]!=' ')
2343                        xx[n++]=xx[i];
2344                    }
2345                    xx[n]='\0';
2346#endif
2347                  }
2348                 
2349                  columnNames = new char * [numberColumns];
2350                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2351                    columnNames[iColumn] =
2352                      strdup(model2->columnName(iColumn).c_str());
2353#ifdef STRIPBLANKS
2354                    char * xx = columnNames[iColumn];
2355                    int i;
2356                    int length = strlen(xx);
2357                    int n=0;
2358                    for (i=0;i<length;i++) {
2359                      if (xx[i]!=' ')
2360                        xx[n++]=xx[i];
2361                    }
2362                    xx[n]='\0';
2363#endif
2364                  }
2365                }
2366                CoinMpsIO writer;
2367                writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
2368                                  model2->getColLower(), model2->getColUpper(),
2369                                  model2->getObjCoefficients(),
2370                                  (const char*) 0 /*integrality*/,
2371                                  model2->getRowLower(), model2->getRowUpper(),
2372                                  columnNames, rowNames);
2373                // Pass in array saying if each variable integer
2374                writer.copyInIntegerInformation(model2->integerInformation());
2375                writer.setObjectiveOffset(model2->objectiveOffset());
2376                writer.writeMps(fileName.c_str(),0,1,1);
2377                if (rowNames) {
2378                  for (iRow=0;iRow<numberRows;iRow++) {
2379                    free(rowNames[iRow]);
2380                  }
2381                  delete [] rowNames;
2382                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2383                    free(columnNames[iColumn]);
2384                  }
2385                  delete [] columnNames;
2386                }
2387#else
2388                model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
2389#endif
2390                if (deleteModel2)
2391                  delete model2;
2392                time2 = CoinCpuTime();
2393                totalTime += time2-time1;
2394                time1=time2;
2395              }
2396            } else {
2397              std::cout<<"** Current model not valid"<<std::endl;
2398            }
2399            break;
2400          case BASISIN:
2401            if (goodModel) {
2402              // get next field
2403              field = CoinReadGetString(argc,argv);
2404              if (field=="$") {
2405                field = parameters[iParam].stringValue();
2406              } else if (field=="EOL") {
2407                parameters[iParam].printString();
2408                break;
2409              } else {
2410                parameters[iParam].setStringValue(field);
2411              }
2412              std::string fileName;
2413              bool canOpen=false;
2414              if (field=="-") {
2415                // stdin
2416                canOpen=true;
2417                fileName = "-";
2418              } else {
2419                if (field[0]=='/'||field[0]=='\\') {
2420                  fileName = field;
2421                } else if (field[0]=='~') {
2422                  char * environVar = getenv("HOME");
2423                  if (environVar) {
2424                    std::string home(environVar);
2425                    field=field.erase(0,1);
2426                    fileName = home+field;
2427                  } else {
2428                    fileName=field;
2429                  }
2430                } else {
2431                  fileName = directory+field;
2432                }
2433                FILE *fp=fopen(fileName.c_str(),"r");
2434                if (fp) {
2435                  // can open - lets go for it
2436                  fclose(fp);
2437                  canOpen=true;
2438                } else {
2439                  std::cout<<"Unable to open file "<<fileName<<std::endl;
2440                }
2441              }
2442              if (canOpen) {
2443                int values = lpSolver->readBasis(fileName.c_str());
2444                if (values==0)
2445                  basisHasValues=-1;
2446                else
2447                  basisHasValues=1;
2448              }
2449            } else {
2450              std::cout<<"** Current model not valid"<<std::endl;
2451            }
2452            break;
2453          case PRIORITYIN:
2454            if (goodModel) {
2455              // get next field
2456              field = CoinReadGetString(argc,argv);
2457              if (field=="$") {
2458                field = parameters[iParam].stringValue();
2459              } else if (field=="EOL") {
2460                parameters[iParam].printString();
2461                break;
2462              } else {
2463                parameters[iParam].setStringValue(field);
2464              }
2465              std::string fileName;
2466              if (field[0]=='/'||field[0]=='\\') {
2467                fileName = field;
2468              } else if (field[0]=='~') {
2469                char * environVar = getenv("HOME");
2470                if (environVar) {
2471                  std::string home(environVar);
2472                  field=field.erase(0,1);
2473                  fileName = home+field;
2474                } else {
2475                  fileName=field;
2476                }
2477              } else {
2478                fileName = directory+field;
2479              }
2480              FILE *fp=fopen(fileName.c_str(),"r");
2481              if (fp) {
2482                // can open - lets go for it
2483                std::string headings[]={"name","number","direction","priority","up","down"};
2484                int order[]={-1,-1,-1,-1,-1,-1};
2485                int nAcross=0;
2486                char line[1000];
2487                int numberColumns = lpSolver->numberColumns();
2488                if (!fgets(line,1000,fp)) {
2489                  std::cout<<"Odd file "<<fileName<<std::endl;
2490                } else {
2491                  char * pos = line;
2492                  char * put = line;
2493                  while (*pos>=' '&&*pos!='\n') {
2494                    if (*pos!=' '&&*pos!='\t') {
2495                      *put=tolower(*pos);
2496                      put++;
2497                    }
2498                  }
2499                  *put='\0';
2500                  pos=line;
2501                  int i;
2502                  bool good=true;
2503                  while (pos) {
2504                    char * comma = strchr(pos,',');
2505                    if (comma)
2506                      *comma='\0';
2507                    for (i=0;i<(int) (sizeof(order)/sizeof(int));i++) {
2508                      if (headings[i]==pos) {
2509                        if (order[i]<0) {
2510                          order[i]=nAcross++;
2511                        } else {
2512                          // duplicate
2513                          good=false;
2514                        }
2515                        break;
2516                      }
2517                    }
2518                    if (i==(int) (sizeof(order)/sizeof(int)))
2519                      good=false;
2520                    if (comma) {
2521                      *comma=',';
2522                      pos=comma+1;
2523                    } else {
2524                      break;
2525                    }
2526                  }
2527                  if (order[0]<0&&order[1]<0)
2528                    good=false;
2529                  if (order[0]>=0&&order[1]>=0)
2530                    good=false;
2531                  if (order[0]>=0&&!lpSolver->lengthNames())
2532                    good=false;
2533                  if (good) {
2534                    char ** columnNames = columnNames = new char * [numberColumns];
2535                    pseudoDown= (double *) malloc(numberColumns*sizeof(double));
2536                    pseudoUp = (double *) malloc(numberColumns*sizeof(double));
2537                    branchDirection = (int *) malloc(numberColumns*sizeof(int));
2538                    priorities= (int *) malloc(numberColumns*sizeof(int));
2539                    int iColumn;
2540                    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2541                      columnNames[iColumn] = 
2542                        strdup(lpSolver->columnName(iColumn).c_str());
2543                      pseudoDown[iColumn]=0.0;
2544                      pseudoUp[iColumn]=0.0;
2545                      branchDirection[iColumn]=0;
2546                      priorities[iColumn]=0;
2547                    }
2548                    int nBadPseudo=0;
2549                    int nBadDir=0;
2550                    int nBadPri=0;
2551                    int nBadName=0;
2552                    int nBadLine=0;
2553                    int nLine=0;
2554                    while (fgets(line,1000,fp)) {
2555                      nLine++;
2556                      iColumn = -1;
2557                      double up =0.0;
2558                      double down=0.0;
2559                      int pri=0;
2560                      int dir=0;
2561                      char * pos = line;
2562                      char * put = line;
2563                      while (*pos>=' '&&*pos!='\n') {
2564                        if (*pos!=' '&&*pos!='\t') {
2565                          *put=tolower(*pos);
2566                          put++;
2567                        }
2568                      }
2569                      *put='\0';
2570                      pos=line;
2571                      for (int i=0;i<nAcross;i++) {
2572                        char * comma = strchr(pos,',');
2573                        if (comma) {
2574                          *comma='\0';
2575                        } else if (i<nAcross-1) {
2576                          nBadLine++;
2577                          break;
2578                        }
2579                        switch (order[i]) {
2580                          // name
2581                        case 0:
2582                          for (iColumn=0;iColumn<numberColumns;iColumn++) {
2583                            if (!strcmp(columnNames[iColumn],pos))
2584                              break;
2585                          }
2586                          if (iColumn==numberColumns)
2587                            iColumn=-1;
2588                          break;
2589                          // number
2590                        case 1:
2591                          iColumn = atoi(pos);
2592                          if (iColumn<0||iColumn>=numberColumns)
2593                            iColumn=-1;
2594                          break;
2595                          // direction
2596                        case 2:
2597                          if (*pos=='D')
2598                            dir=-1;
2599                          else if (*pos=='U')
2600                            dir=1;
2601                          else if (*pos=='N')
2602                            dir=0;
2603                          else if (*pos=='1'&&*(pos+1)=='\0')
2604                            dir=1;
2605                          else if (*pos=='0'&&*(pos+1)=='\0')
2606                            dir=0;
2607                          else if (*pos=='1'&&*(pos+1)=='1'&&*(pos+2)=='\0')
2608                            dir=-1;
2609                          else
2610                            dir=-2; // bad
2611                          break;
2612                          // priority
2613                        case 3:
2614                          pri=atoi(pos);
2615                          break;
2616                          // up
2617                        case 4:
2618                          up = atof(pos);
2619                          break;
2620                          // down
2621                        case 5:
2622                          down = atof(pos);
2623                          break;
2624                        }
2625                        if (comma) {
2626                          *comma=',';
2627                          pos=comma+1;
2628                        }
2629                      }
2630                      if (iColumn>=0) {
2631                        if (down<0.0) {
2632                          nBadPseudo++;
2633                          down=0.0;
2634                        }
2635                        if (up<0.0) {
2636                          nBadPseudo++;
2637                          up=0.0;
2638                        }
2639                        if (!up)
2640                          up=down;
2641                        if (!down)
2642                          down=up;
2643                        if (dir<-1||dir>1) {
2644                          nBadDir++;
2645                          dir=0;
2646                        }
2647                        if (pri<0) {
2648                          nBadPri++;
2649                          pri=0;
2650                        }
2651                        pseudoDown[iColumn]=down;
2652                        pseudoUp[iColumn]=up;
2653                        branchDirection[iColumn]=dir;
2654                        priorities[iColumn]=pri;
2655                      } else {
2656                        nBadName++;
2657                      }
2658                    }
2659                    if (!noPrinting) {
2660                      printf("%d fields and %d records",nAcross,nLine);
2661                      if (nBadPseudo)
2662                        printf(" %d bad pseudo costs",nBadPseudo);
2663                      if (nBadDir)
2664                        printf(" %d bad directions",nBadDir);
2665                      if (nBadPri)
2666                        printf(" %d bad priorities",nBadPri);
2667                      if (nBadName)
2668                        printf(" ** %d records did not match on name/sequence",nBadName);
2669                      printf("\n");
2670                    }
2671                    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2672                      free(columnNames[iColumn]);
2673                    }
2674                    delete [] columnNames;
2675                  } else {
2676                    std::cout<<"Duplicate or unknown keyword - or name/number fields wrong"<<line<<std::endl;
2677                  }
2678                }
2679                fclose(fp);
2680              } else {
2681                std::cout<<"Unable to open file "<<fileName<<std::endl;
2682              }
2683            } else {
2684              std::cout<<"** Current model not valid"<<std::endl;
2685            }
2686            break;
2687          case DEBUG:
2688            if (goodModel) {
2689              delete [] debugValues;
2690              debugValues=NULL;
2691              // get next field
2692              field = CoinReadGetString(argc,argv);
2693              if (field=="$") {
2694                field = parameters[iParam].stringValue();
2695              } else if (field=="EOL") {
2696                parameters[iParam].printString();
2697                break;
2698              } else {
2699                parameters[iParam].setStringValue(field);
2700                debugFile=field;
2701                if (debugFile=="create"||
2702                    debugFile=="createAfterPre") {
2703                  printf("Will create a debug file so this run should be a good one\n");
2704                  break;
2705                }
2706              }
2707              std::string fileName;
2708              if (field[0]=='/'||field[0]=='\\') {
2709                fileName = field;
2710              } else if (field[0]=='~') {
2711                char * environVar = getenv("HOME");
2712                if (environVar) {
2713                  std::string home(environVar);
2714                  field=field.erase(0,1);
2715                  fileName = home+field;
2716                } else {
2717                  fileName=field;
2718                }
2719              } else {
2720                fileName = directory+field;
2721              }
2722              FILE *fp=fopen(fileName.c_str(),"rb");
2723              if (fp) {
2724                // can open - lets go for it
2725                int numRows;
2726                double obj;
2727                fread(&numRows,sizeof(int),1,fp);
2728                fread(&numberDebugValues,sizeof(int),1,fp);
2729                fread(&obj,sizeof(double),1,fp);
2730                debugValues = new double[numberDebugValues+numRows];
2731                fread(debugValues,sizeof(double),numRows,fp);
2732                fread(debugValues,sizeof(double),numRows,fp);
2733                fread(debugValues,sizeof(double),numberDebugValues,fp);
2734                printf("%d doubles read into debugValues\n",numberDebugValues);
2735                fclose(fp);
2736              } else {
2737                std::cout<<"Unable to open file "<<fileName<<std::endl;
2738              }
2739            } else {
2740              std::cout<<"** Current model not valid"<<std::endl;
2741            }
2742            break;
2743          case PRINTMASK:
2744            // get next field
2745            {
2746              std::string name = CoinReadGetString(argc,argv);
2747              if (name!="EOL") {
2748                parameters[iParam].setStringValue(name);
2749                printMask = name;
2750              } else {
2751                parameters[iParam].printString();
2752              }
2753            }
2754            break;
2755          case BASISOUT:
2756            if (goodModel) {
2757              // get next field
2758              field = CoinReadGetString(argc,argv);
2759              if (field=="$") {
2760                field = parameters[iParam].stringValue();
2761              } else if (field=="EOL") {
2762                parameters[iParam].printString();
2763                break;
2764              } else {
2765                parameters[iParam].setStringValue(field);
2766              }
2767              std::string fileName;
2768              bool canOpen=false;
2769              if (field[0]=='/'||field[0]=='\\') {
2770                fileName = field;
2771              } else if (field[0]=='~') {
2772                char * environVar = getenv("HOME");
2773                if (environVar) {
2774                  std::string home(environVar);
2775                  field=field.erase(0,1);
2776                  fileName = home+field;
2777                } else {
2778                  fileName=field;
2779                }
2780              } else {
2781                fileName = directory+field;
2782              }
2783              FILE *fp=fopen(fileName.c_str(),"w");
2784              if (fp) {
2785                // can open - lets go for it
2786                fclose(fp);
2787                canOpen=true;
2788              } else {
2789                std::cout<<"Unable to open file "<<fileName<<std::endl;
2790              }
2791              if (canOpen) {
2792                ClpSimplex * model2 = lpSolver;
2793                model2->writeBasis(fileName.c_str(),outputFormat>1,outputFormat-2);
2794                time2 = CoinCpuTime();
2795                totalTime += time2-time1;
2796                time1=time2;
2797              }
2798            } else {
2799              std::cout<<"** Current model not valid"<<std::endl;
2800            }
2801            break;
2802          case SAVE:
2803            {
2804              // get next field
2805              field = CoinReadGetString(argc,argv);
2806              if (field=="$") {
2807                field = parameters[iParam].stringValue();
2808              } else if (field=="EOL") {
2809                parameters[iParam].printString();
2810                break;
2811              } else {
2812                parameters[iParam].setStringValue(field);
2813              }
2814              std::string fileName;
2815              bool canOpen=false;
2816              if (field[0]=='/'||field[0]=='\\') {
2817                fileName = field;
2818              } else if (field[0]=='~') {
2819                char * environVar = getenv("HOME");
2820                if (environVar) {
2821                  std::string home(environVar);
2822                  field=field.erase(0,1);
2823                  fileName = home+field;
2824                } else {
2825                  fileName=field;
2826                }
2827              } else {
2828                fileName = directory+field;
2829              }
2830              FILE *fp=fopen(fileName.c_str(),"wb");
2831              if (fp) {
2832                // can open - lets go for it
2833                fclose(fp);
2834                canOpen=true;
2835              } else {
2836                std::cout<<"Unable to open file "<<fileName<<std::endl;
2837              }
2838              if (canOpen) {
2839                int status;
2840                // If presolve on then save presolved
2841                bool deleteModel2=false;
2842                ClpSimplex * model2 = lpSolver;
2843                if (preSolve) {
2844                  ClpPresolve pinfo;
2845                  double presolveTolerance = 
2846                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
2847                  model2 = 
2848                    pinfo.presolvedModel(*lpSolver,presolveTolerance,
2849                                         false,preSolve);
2850                  if (model2) {
2851                    printf("Saving presolved model on %s\n",
2852                           fileName.c_str());
2853                    deleteModel2=true;
2854                  } else {
2855                    printf("Presolved model looks infeasible - saving original on %s\n",
2856                           fileName.c_str());
2857                    deleteModel2=false;
2858                    model2 = lpSolver;
2859
2860                  }
2861                } else {
2862                  printf("Saving model on %s\n",
2863                           fileName.c_str());
2864                }
2865                status =model2->saveModel(fileName.c_str());
2866                if (deleteModel2)
2867                  delete model2;
2868                if (!status) {
2869                  goodModel=true;
2870                  time2 = CoinCpuTime();
2871                  totalTime += time2-time1;
2872                  time1=time2;
2873                } else {
2874                  // errors
2875                  std::cout<<"There were errors on output"<<std::endl;
2876                }
2877              }
2878            }
2879            break;
2880          case RESTORE:
2881            {
2882              // get next field
2883              field = CoinReadGetString(argc,argv);
2884              if (field=="$") {
2885                field = parameters[iParam].stringValue();
2886              } else if (field=="EOL") {
2887                parameters[iParam].printString();
2888                break;
2889              } else {
2890                parameters[iParam].setStringValue(field);
2891              }
2892              std::string fileName;
2893              bool canOpen=false;
2894              if (field[0]=='/'||field[0]=='\\') {
2895                fileName = field;
2896              } else if (field[0]=='~') {
2897                char * environVar = getenv("HOME");
2898                if (environVar) {
2899                  std::string home(environVar);
2900                  field=field.erase(0,1);
2901                  fileName = home+field;
2902                } else {
2903                  fileName=field;
2904                }
2905              } else {
2906                fileName = directory+field;
2907              }
2908              FILE *fp=fopen(fileName.c_str(),"rb");
2909              if (fp) {
2910                // can open - lets go for it
2911                fclose(fp);
2912                canOpen=true;
2913              } else {
2914                std::cout<<"Unable to open file "<<fileName<<std::endl;
2915              }
2916              if (canOpen) {
2917                int status =lpSolver->restoreModel(fileName.c_str());
2918                if (!status) {
2919                  goodModel=true;
2920                  time2 = CoinCpuTime();
2921                  totalTime += time2-time1;
2922                  time1=time2;
2923                } else {
2924                  // errors
2925                  std::cout<<"There were errors on input"<<std::endl;
2926                }
2927              }
2928            }
2929            break;
2930          case MAXIMIZE:
2931            lpSolver->setOptimizationDirection(-1);
2932            break;
2933          case MINIMIZE:
2934            lpSolver->setOptimizationDirection(1);
2935            break;
2936          case ALLSLACK:
2937            lpSolver->allSlackBasis(true);
2938            break;
2939          case REVERSE:
2940            if (goodModel) {
2941              int iColumn;
2942              int numberColumns=lpSolver->numberColumns();
2943              double * dualColumnSolution = 
2944                lpSolver->dualColumnSolution();
2945              ClpObjective * obj = lpSolver->objectiveAsObject();
2946              assert(dynamic_cast<ClpLinearObjective *> (obj));
2947              double offset;
2948              double * objective = obj->gradient(NULL,NULL,offset,true);
2949              for (iColumn=0;iColumn<numberColumns;iColumn++) {
2950                dualColumnSolution[iColumn] = dualColumnSolution[iColumn];
2951                objective[iColumn] = -objective[iColumn];
2952              }
2953              int iRow;
2954              int numberRows=lpSolver->numberRows();
2955              double * dualRowSolution = 
2956                lpSolver->dualRowSolution();
2957              for (iRow=0;iRow<numberRows;iRow++) 
2958                dualRowSolution[iRow] = dualRowSolution[iRow];
2959            }
2960            break;
2961          case DIRECTORY:
2962            {
2963              std::string name = CoinReadGetString(argc,argv);
2964              if (name!="EOL") {
2965                int length=name.length();
2966                if (name[length-1]=='/'||name[length-1]=='\\')
2967                  directory=name;
2968                else
2969                  directory = name+"/";
2970                parameters[iParam].setStringValue(directory);
2971              } else {
2972                parameters[iParam].printString();
2973              }
2974            }
2975            break;
2976          case STDIN:
2977            CbcOrClpRead_mode=-1;
2978            break;
2979          case NETLIB_DUAL:
2980          case NETLIB_EITHER:
2981          case NETLIB_BARRIER:
2982          case NETLIB_PRIMAL:
2983          case NETLIB_TUNE:
2984            {
2985              // create fields for unitTest
2986              const char * fields[4];
2987              int nFields=2;
2988              fields[0]="fake main from unitTest";
2989              fields[1]="-netlib";
2990              if (directory!=defaultDirectory) {
2991                fields[2]="-netlibDir";
2992                fields[3]=directory.c_str();
2993                nFields=4;
2994              }
2995              int algorithm;
2996              if (type==NETLIB_DUAL) {
2997                std::cerr<<"Doing netlib with dual algorithm"<<std::endl;
2998                algorithm =0;
2999              } else if (type==NETLIB_BARRIER) {
3000                std::cerr<<"Doing netlib with barrier algorithm"<<std::endl;
3001                algorithm =2;
3002              } else if (type==NETLIB_EITHER) {
3003                std::cerr<<"Doing netlib with dual or primal algorithm"<<std::endl;
3004                algorithm =3;
3005              } else if (type==NETLIB_TUNE) {
3006                std::cerr<<"Doing netlib with best algorithm!"<<std::endl;
3007                algorithm =5;
3008                // uncomment next to get active tuning
3009                // algorithm=6;
3010              } else {
3011                std::cerr<<"Doing netlib with primal agorithm"<<std::endl;
3012                algorithm=1;
3013              }
3014              int specialOptions = lpSolver->specialOptions();
3015              lpSolver->setSpecialOptions(0);
3016              mainTest(nFields,fields,algorithm,*lpSolver,
3017                       (preSolve!=0),specialOptions);
3018            }
3019            break;
3020          case UNITTEST:
3021            {
3022              // create fields for unitTest
3023              const char * fields[3];
3024              int nFields=1;
3025              fields[0]="fake main from unitTest";
3026              if (directory!=defaultDirectory) {
3027                fields[1]="-mpsDir";
3028                fields[2]=directory.c_str();
3029                nFields=3;
3030              }
3031              mainTest(nFields,fields,false,*lpSolver,(preSolve!=0),
3032                       false);
3033            }
3034            break;
3035          case FAKEBOUND:
3036            if (goodModel) {
3037              // get bound
3038              double value = CoinReadGetDoubleField(argc,argv,&valid);
3039              if (!valid) {
3040                std::cout<<"Setting "<<parameters[iParam].name()<<
3041                  " to DEBUG "<<value<<std::endl;
3042                int iRow;
3043                int numberRows=lpSolver->numberRows();
3044                double * rowLower = lpSolver->rowLower();
3045                double * rowUpper = lpSolver->rowUpper();
3046                for (iRow=0;iRow<numberRows;iRow++) {
3047                  // leave free ones for now
3048                  if (rowLower[iRow]>-1.0e20||rowUpper[iRow]<1.0e20) {
3049                    rowLower[iRow]=CoinMax(rowLower[iRow],-value);
3050                    rowUpper[iRow]=CoinMin(rowUpper[iRow],value);
3051                  }
3052                }
3053                int iColumn;
3054                int numberColumns=lpSolver->numberColumns();
3055                double * columnLower = lpSolver->columnLower();
3056                double * columnUpper = lpSolver->columnUpper();
3057                for (iColumn=0;iColumn<numberColumns;iColumn++) {
3058                  // leave free ones for now
3059                  if (columnLower[iColumn]>-1.0e20||
3060                      columnUpper[iColumn]<1.0e20) {
3061                    columnLower[iColumn]=CoinMax(columnLower[iColumn],-value);
3062                    columnUpper[iColumn]=CoinMin(columnUpper[iColumn],value);
3063                  }
3064                }
3065              } else if (valid==1) {
3066                abort();
3067              } else {
3068                std::cout<<"enter value for "<<parameters[iParam].name()<<
3069                  std::endl;
3070              }
3071            }
3072            break;
3073          case REALLY_SCALE:
3074            if (goodModel) {
3075              ClpSimplex newModel(*lpSolver,
3076                                  lpSolver->scalingFlag());
3077              printf("model really really scaled\n");
3078              *lpSolver=newModel;
3079            }
3080            break;
3081          case HELP:
3082            std::cout<<"Coin Solver version "<<CBCVERSION
3083                     <<", build "<<__DATE__<<std::endl;
3084            std::cout<<"Non default values:-"<<std::endl;
3085            std::cout<<"Perturbation "<<lpSolver->perturbation()<<" (default 100)"
3086                     <<std::endl;
3087            CoinReadPrintit(
3088                    "Presolve being done with 5 passes\n\
3089Dual steepest edge steep/partial on matrix shape and factorization density\n\
3090Clpnnnn taken out of messages\n\
3091If Factorization frequency default then done on size of matrix\n\n\
3092(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
3093You can switch to interactive mode at any time so\n\
3094clp watson.mps -scaling off -primalsimplex\nis the same as\n\
3095clp watson.mps -\nscaling off\nprimalsimplex"
3096                    );
3097            break;
3098          case SOLUTION:
3099            if (goodModel) {
3100              // get next field
3101              field = CoinReadGetString(argc,argv);
3102              if (field=="$") {
3103                field = parameters[iParam].stringValue();
3104              } else if (field=="EOL") {
3105                parameters[iParam].printString();
3106                break;
3107              } else {
3108                parameters[iParam].setStringValue(field);
3109              }
3110              std::string fileName;
3111              FILE *fp=NULL;
3112              if (field=="-"||field=="EOL"||field=="stdout") {
3113                // stdout
3114                fp=stdout;
3115              } else if (field=="stderr") {
3116                // stderr
3117                fp=stderr;
3118              } else {
3119                if (field[0]=='/'||field[0]=='\\') {
3120                  fileName = field;
3121                } else if (field[0]=='~') {
3122                  char * environVar = getenv("HOME");
3123                  if (environVar) {
3124                    std::string home(environVar);
3125                    field=field.erase(0,1);
3126                    fileName = home+field;
3127                  } else {
3128                    fileName=field;
3129                  }
3130                } else {
3131                  fileName = directory+field;
3132                }
3133                fp=fopen(fileName.c_str(),"w");
3134              }
3135              if (fp) {
3136                // make fancy later on
3137                int iRow;
3138                int numberRows=lpSolver->numberRows();
3139                double * dualRowSolution = lpSolver->dualRowSolution();
3140                double * primalRowSolution = 
3141                  lpSolver->primalRowSolution();
3142                double * rowLower = lpSolver->rowLower();
3143                double * rowUpper = lpSolver->rowUpper();
3144                double primalTolerance = lpSolver->primalTolerance();
3145                char format[6];
3146                sprintf(format,"%%-%ds",CoinMax(lengthName,8));
3147                bool doMask = (printMask!=""&&lengthName);
3148                if (printMode>2) {
3149                  for (iRow=0;iRow<numberRows;iRow++) {
3150                    int type=printMode-3;
3151                    if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
3152                        primalRowSolution[iRow]<rowLower[iRow]-primalTolerance) {
3153                      fprintf(fp,"** ");
3154                      type=2;
3155                    } else if (fabs(primalRowSolution[iRow])>1.0e-8) {
3156                      type=1;
3157                    } else if (numberRows<50) {
3158                      type=3;
3159                    }
3160                    if (doMask&&!maskMatches(printMask,rowNames[iRow]))
3161                      type =0;
3162                    if (type) {
3163                      fprintf(fp,"%7d ",iRow);
3164                      if (lengthName)
3165                        fprintf(fp,format,rowNames[iRow].c_str());
3166                      fprintf(fp,"%15.8g        %15.8g\n",primalRowSolution[iRow],
3167                              dualRowSolution[iRow]);
3168                    }
3169                  }
3170                }
3171                int iColumn;
3172                int numberColumns=lpSolver->numberColumns();
3173                double * dualColumnSolution = 
3174                  lpSolver->dualColumnSolution();
3175                double * primalColumnSolution = 
3176                  lpSolver->primalColumnSolution();
3177                double * columnLower = lpSolver->columnLower();
3178                double * columnUpper = lpSolver->columnUpper();
3179                if (printMode!=2) {
3180                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3181                    int type=(printMode>3) ? 1 :0;
3182                    if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
3183                        primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance) {
3184                      fprintf(fp,"** ");
3185                      type=2;
3186                    } else if (fabs(primalColumnSolution[iColumn])>1.0e-8) {
3187                      type=1;
3188                    } else if (numberColumns<50) {
3189                      type=3;
3190                    }
3191                    // see if integer
3192                    if ((!lpSolver->isInteger(iColumn)||fabs(primalColumnSolution[iColumn])<1.0e-8)
3193                         &&printMode==1)
3194                      type=0;
3195                    if (doMask&&!maskMatches(printMask,columnNames[iColumn]))
3196                      type =0;
3197                    if (type) {
3198                      fprintf(fp,"%7d ",iColumn);
3199                      if (lengthName)
3200                        fprintf(fp,format,columnNames[iColumn].c_str());
3201                      fprintf(fp,"%15.8g        %15.8g\n",
3202                              primalColumnSolution[iColumn],
3203                              dualColumnSolution[iColumn]);
3204                    }
3205                  }
3206                } else {
3207                  // special format suitable for OsiRowCutDebugger
3208                  int n=0;
3209                  bool comma=false;
3210                  bool newLine=false;
3211                  fprintf(fp,"\tint intIndicesV[]={\n");
3212                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3213                    if(primalColumnSolution[iColumn]>0.5&&model.solver()->isInteger(iColumn)) {
3214                      if (comma)
3215                        fprintf(fp,",");
3216                      if (newLine)
3217                        fprintf(fp,"\n");
3218                      fprintf(fp,"%d ",iColumn);
3219                      comma=true;
3220                      newLine=false;
3221                      n++;
3222                      if (n==10) {
3223                        n=0;
3224                        newLine=true;
3225                      }
3226                    }
3227                  }
3228                  fprintf(fp,"};\n");
3229                  n=0;
3230                  comma=false;
3231                  newLine=false;
3232                  fprintf(fp,"\tdouble intSolnV[]={\n");
3233                  for ( iColumn=0;iColumn<numberColumns;iColumn++) {
3234                    if(primalColumnSolution[iColumn]>0.5&&model.solver()->isInteger(iColumn)) {
3235                      if (comma)
3236                        fprintf(fp,",");
3237                      if (newLine)
3238                        fprintf(fp,"\n");
3239                      int value = (int) (primalColumnSolution[iColumn]+0.5);
3240                      fprintf(fp,"%d. ",value);
3241                      comma=true;
3242                      newLine=false;
3243                      n++;
3244                      if (n==10) {
3245                        n=0;
3246                        newLine=true;
3247                      }
3248                    }
3249                  }
3250                  fprintf(fp,"};\n");
3251                }
3252                if (fp!=stdout)
3253                  fclose(fp);
3254              } else {
3255                std::cout<<"Unable to open file "<<fileName<<std::endl;
3256              }
3257            } else {
3258              std::cout<<"** Current model not valid"<<std::endl;
3259             
3260            }
3261            break;
3262          case SAVESOL:
3263            if (goodModel) {
3264              // get next field
3265              field = CoinReadGetString(argc,argv);
3266              if (field=="$") {
3267                field = parameters[iParam].stringValue();
3268              } else if (field=="EOL") {
3269                parameters[iParam].printString();
3270                break;
3271              } else {
3272                parameters[iParam].setStringValue(field);
3273              }
3274              std::string fileName;
3275              if (field[0]=='/'||field[0]=='\\') {
3276                fileName = field;
3277              } else if (field[0]=='~') {
3278                char * environVar = getenv("HOME");
3279                if (environVar) {
3280                  std::string home(environVar);
3281                  field=field.erase(0,1);
3282                  fileName = home+field;
3283                } else {
3284                  fileName=field;
3285                }
3286              } else {
3287                fileName = directory+field;
3288              }
3289              saveSolution(lpSolver,fileName);
3290            } else {
3291              std::cout<<"** Current model not valid"<<std::endl;
3292             
3293            }
3294            break;
3295          case DUMMY:
3296            break;
3297          default:
3298            abort();
3299          }
3300        } 
3301      } else if (!numberMatches) {
3302        std::cout<<"No match for "<<field<<" - ? for list of commands"
3303                 <<std::endl;
3304      } else if (numberMatches==1) {
3305        if (!numberQuery) {
3306          std::cout<<"Short match for "<<field<<" - completion: ";
3307          std::cout<<parameters[firstMatch].matchName()<<std::endl;
3308        } else if (numberQuery) {
3309          std::cout<<parameters[firstMatch].matchName()<<" : ";
3310          std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
3311          if (numberQuery>=2) 
3312            parameters[firstMatch].printLongHelp();
3313        }
3314      } else {
3315        if (!numberQuery) 
3316          std::cout<<"Multiple matches for "<<field<<" - possible completions:"
3317                   <<std::endl;
3318        else
3319          std::cout<<"Completions of "<<field<<":"<<std::endl;
3320        for ( iParam=0; iParam<numberParameters; iParam++ ) {
3321          int match = parameters[iParam].matches(field);
3322          if (match&&parameters[iParam].displayThis()) {
3323            std::cout<<parameters[iParam].matchName();
3324            if (numberQuery>=2) 
3325              std::cout<<" : "<<parameters[iParam].shortHelp();
3326            std::cout<<std::endl;
3327          }
3328        }
3329      }
3330    }
3331  }
3332  // By now all memory should be freed
3333#ifdef DMALLOC
3334  dmalloc_log_unfreed();
3335  dmalloc_shutdown();
3336#endif
3337  return 0;
3338}   
3339static void breakdown(const char * name, int numberLook, const double * region)
3340{
3341  double range[] = {
3342    -COIN_DBL_MAX,
3343    -1.0e15,-1.0e11,-1.0e8,-1.0e5,-1.0e4,-1.0e3,-1.0e2,-1.0e1,
3344    -1.0,
3345    -1.0e-1,-1.0e-2,-1.0e-3,-1.0e-4,-1.0e-5,-1.0e-8,-1.0e-11,-1.0e-15,
3346    0.0,
3347    1.0e-15,1.0e-11,1.0e-8,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,
3348    1.0,
3349    1.0e1,1.0e2,1.0e3,1.0e4,1.0e5,1.0e8,1.0e11,1.0e15,
3350    COIN_DBL_MAX};
3351  int nRanges = (int) (sizeof(range)/sizeof(double));
3352  int * number = new int[nRanges];
3353  memset(number,0,nRanges*sizeof(int));
3354  int * numberExact = new int[nRanges];
3355  memset(numberExact,0,nRanges*sizeof(int));
3356  int i;
3357  for ( i=0;i<numberLook;i++) {
3358    double value = region[i];
3359    for (int j=0;j<nRanges;j++) {
3360      if (value==range[j]) {
3361        numberExact[j]++;
3362        break;
3363      } else if (value<range[j]) {
3364        number[j]++;
3365        break;
3366      }
3367    }
3368  }
3369  printf("\n%s has %d entries\n",name,numberLook);
3370  for (i=0;i<nRanges;i++) {
3371    if (number[i]) 
3372      printf("%d between %g and %g",number[i],range[i-1],range[i]);
3373    if (numberExact[i]) {
3374      if (number[i])
3375        printf(", ");
3376      printf("%d exactly at %g",numberExact[i],range[i]);
3377    }
3378    if (number[i]+numberExact[i])
3379      printf("\n");
3380  }
3381  delete [] number;
3382  delete [] numberExact;
3383}
3384static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
3385{
3386  int numberColumns = originalModel->numberColumns();
3387  const char * integerInformation  = originalModel->integerInformation(); 
3388  const double * columnLower = originalModel->columnLower();
3389  const double * columnUpper = originalModel->columnUpper();
3390  int numberIntegers=0;
3391  int numberBinary=0;
3392  int iRow,iColumn;
3393  if (integerInformation) {
3394    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3395      if (integerInformation[iColumn]) {
3396        if (columnUpper[iColumn]>columnLower[iColumn]) {
3397          numberIntegers++;
3398          if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==1) 
3399            numberBinary++;
3400        }
3401      }
3402    }
3403  }
3404  numberColumns = model->numberColumns();
3405  int numberRows = model->numberRows();
3406  columnLower = model->columnLower();
3407  columnUpper = model->columnUpper();
3408  const double * rowLower = model->rowLower();
3409  const double * rowUpper = model->rowUpper();
3410  const double * objective = model->objective();
3411  CoinPackedMatrix * matrix = model->matrix();
3412  CoinBigIndex numberElements = matrix->getNumElements();
3413  const int * columnLength = matrix->getVectorLengths();
3414  //const CoinBigIndex * columnStart = matrix->getVectorStarts();
3415  const double * elementByColumn = matrix->getElements();
3416  int * number = new int[numberRows+1];
3417  memset(number,0,(numberRows+1)*sizeof(int));
3418  int numberObjSingletons=0;
3419  /* cType
3420     0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
3421     8 0/1
3422  */ 
3423  int cType[9];
3424  std::string cName[]={"0.0->inf,","0.0->up,","lo->inf,","lo->up,","free,","fixed,","-inf->0.0,",
3425                       "-inf->up,","0.0->1.0"};
3426  int nObjective=0;
3427  memset(cType,0,sizeof(cType));
3428  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3429    int length=columnLength[iColumn];
3430    if (length==1&&objective[iColumn])
3431      numberObjSingletons++;
3432    number[length]++;
3433    if (objective[iColumn])
3434      nObjective++;
3435    if (columnLower[iColumn]>-1.0e20) {
3436      if (columnLower[iColumn]==0.0) {
3437        if (columnUpper[iColumn]>1.0e20)
3438          cType[0]++;
3439        else if (columnUpper[iColumn]==1.0)
3440          cType[8]++;
3441        else if (columnUpper[iColumn]==0.0)
3442          cType[5]++;
3443        else
3444          cType[1]++;
3445      } else {
3446        if (columnUpper[iColumn]>1.0e20) 
3447          cType[2]++;
3448        else if (columnUpper[iColumn]==columnLower[iColumn])
3449          cType[5]++;
3450        else
3451          cType[3]++;
3452      }
3453    } else {
3454      if (columnUpper[iColumn]>1.0e20) 
3455        cType[4]++;
3456      else if (columnUpper[iColumn]==0.0) 
3457        cType[6]++;
3458      else
3459        cType[7]++;
3460    }
3461  }
3462  /* rType
3463     0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
3464     7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
3465  */ 
3466  int rType[13];
3467  std::string rName[]={"E 0.0,","E 1.0,","E -1.0,","E other,","G 0.0,","G 1.0,","G other,",
3468                       "L 0.0,","L 1.0,","L other,","Range 0.0->1.0,","Range other,","Free"};
3469  memset(rType,0,sizeof(rType));
3470  for (iRow=0;iRow<numberRows;iRow++) {
3471    if (rowLower[iRow]>-1.0e20) {
3472      if (rowLower[iRow]==0.0) {
3473        if (rowUpper[iRow]>1.0e20)
3474          rType[4]++;
3475        else if (rowUpper[iRow]==1.0)
3476          rType[10]++;
3477        else if (rowUpper[iRow]==0.0)
3478          rType[0]++;
3479        else
3480          rType[11]++;
3481      } else if (rowLower[iRow]==1.0) {
3482        if (rowUpper[iRow]>1.0e20) 
3483          rType[5]++;
3484        else if (rowUpper[iRow]==rowLower[iRow])
3485          rType[1]++;
3486        else
3487          rType[11]++;
3488      } else if (rowLower[iRow]==-1.0) {
3489        if (rowUpper[iRow]>1.0e20) 
3490          rType[6]++;
3491        else if (rowUpper[iRow]==rowLower[iRow])
3492          rType[2]++;
3493        else
3494          rType[11]++;
3495      } else {
3496        if (rowUpper[iRow]>1.0e20) 
3497          rType[6]++;
3498        else if (rowUpper[iRow]==rowLower[iRow])
3499          rType[3]++;
3500        else
3501          rType[11]++;
3502      }
3503    } else {
3504      if (rowUpper[iRow]>1.0e20) 
3505        rType[12]++;
3506      else if (rowUpper[iRow]==0.0) 
3507        rType[7]++;
3508      else if (rowUpper[iRow]==1.0) 
3509        rType[8]++;
3510      else
3511        rType[9]++;
3512    }
3513  }
3514  // Basic statistics
3515  printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
3516         numberRows,numberColumns,nObjective,numberElements);
3517  if (number[0]+number[1]) {
3518    printf("There are ");
3519    if (numberObjSingletons)
3520      printf("%d singletons with objective ",numberObjSingletons);
3521    int numberNoObj = number[1]-numberObjSingletons;
3522    if (numberNoObj)
3523      printf("%d singletons with no objective ",numberNoObj);
3524    if (number[0])
3525      printf("** %d columns have no entries",number[0]);
3526    printf("\n");
3527  }
3528  printf("Column breakdown:\n");
3529  int k;
3530  for (k=0;k<(int) (sizeof(cType)/sizeof(int));k++) {
3531    printf("%d of type %s ",cType[k],cName[k].c_str());
3532    if (((k+1)%3)==0)
3533      printf("\n");
3534  }
3535  if ((k%3)!=0)
3536    printf("\n");
3537  printf("Row breakdown:\n");
3538  for (k=0;k<(int) (sizeof(rType)/sizeof(int));k++) {
3539    printf("%d of type %s ",rType[k],rName[k].c_str());
3540    if (((k+1)%3)==0)
3541      printf("\n");
3542  }
3543  if ((k%3)!=0)
3544    printf("\n");
3545  if (model->logLevel()<2)
3546    return ;
3547  int kMax = model->logLevel()>3 ? 1000000 : 10;
3548  k=0;
3549  for (iRow=1;iRow<=numberRows;iRow++) {
3550    if (number[iRow]) {
3551      k++;
3552      printf("%d columns have %d entries\n",number[iRow],iRow);
3553      if (k==kMax)
3554        break;
3555    }
3556  }
3557  if (k<numberRows) {
3558    int kk=k;
3559    k=0;
3560    for (iRow=numberRows;iRow>=1;iRow--) {
3561      if (number[iRow]) {
3562        k++;
3563        if (k==kMax)
3564          break;
3565      }
3566    }
3567    if (k>kk) {
3568      printf("\n    .........\n\n");
3569      iRow=k;
3570      k=0;
3571      for (;iRow<numberRows;iRow++) {
3572        if (number[iRow]) {
3573          k++;
3574          printf("%d columns have %d entries\n",number[iRow],iRow);
3575          if (k==kMax)
3576            break;
3577        }
3578      }
3579    }
3580  }
3581  delete [] number;
3582  printf("\n\n");
3583  // get row copy
3584  CoinPackedMatrix rowCopy = *matrix;
3585  rowCopy.reverseOrdering();
3586  //const int * column = rowCopy.getIndices();
3587  const int * rowLength = rowCopy.getVectorLengths();
3588  //const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
3589  //const double * element = rowCopy.getElements();
3590  number = new int[numberColumns+1];
3591  memset(number,0,(numberColumns+1)*sizeof(int));
3592  for (iRow=0;iRow<numberRows;iRow++) {
3593    int length=rowLength[iRow];
3594    number[length]++;
3595  }
3596  if (number[0])
3597    printf("** %d rows have no entries\n",number[0]);
3598  k=0;
3599  for (iColumn=1;iColumn<=numberColumns;iColumn++) {
3600    if (number[iColumn]) {
3601      k++;
3602      printf("%d rows have %d entries\n",number[iColumn],iColumn);
3603      if (k==kMax)
3604        break;
3605    }
3606  }
3607  if (k<numberColumns) {
3608    int kk=k;
3609    k=0;
3610    for (iColumn=numberColumns;iColumn>=1;iColumn--) {
3611      if (number[iColumn]) {
3612        k++;
3613        if (k==kMax)
3614          break;
3615      }
3616    }
3617    if (k>kk) {
3618      printf("\n    .........\n\n");
3619      iColumn=k;
3620      k=0;
3621      for (;iColumn<numberColumns;iColumn++) {
3622        if (number[iColumn]) {
3623          k++;
3624          printf("%d rows have %d entries\n",number[iColumn],iColumn);
3625          if (k==kMax)
3626            break;
3627        }
3628      }
3629    }
3630  }
3631  delete [] number;
3632  // Now do breakdown of ranges
3633  breakdown("Elements",numberElements,elementByColumn);
3634  breakdown("RowLower",numberRows,rowLower);
3635  breakdown("RowUpper",numberRows,rowUpper);
3636  breakdown("ColumnLower",numberColumns,columnLower);
3637  breakdown("ColumnUpper",numberColumns,columnUpper);
3638  breakdown("Objective",numberColumns,objective);
3639}
3640static bool maskMatches(std::string & mask, std::string & check)
3641{
3642  // back to char as I am old fashioned
3643  const char * maskC = mask.c_str();
3644  const char * checkC = check.c_str();
3645  int length = strlen(maskC);
3646  int lengthCheck;
3647  for (lengthCheck=length-1;lengthCheck>=0;lengthCheck--) {
3648    if (maskC[lengthCheck]!='*')
3649      break;
3650  }
3651  lengthCheck++;
3652  int lengthC = strlen(checkC);
3653  if (lengthC>length)
3654    return false; // can't be true
3655  if (lengthC<lengthCheck) {
3656    // last lot must be blank for match
3657    for (int i=lengthC;i<lengthCheck;i++) {
3658      if (maskC[i]!=' ')
3659        return false;
3660    }
3661  }
3662  // need only check this much
3663  lengthC = CoinMin(lengthC,lengthCheck);
3664  for (int i=0;i<lengthC;i++) {
3665    if (maskC[i]!='*'&&maskC[i]!=checkC[i])
3666      return false;
3667  }
3668  return true; // matches
3669}
3670/*
3671  Version 1.00.00 November 16 2005.
3672  This is to stop me (JJF) messing about too much.
3673  Tuning changes should be noted here.
3674  The testing next version may be activated by CBC_NEXT_VERSION
3675  This applies to OsiClp, Clp etc
3676  Version 1.00.01 November 24 2005
3677  Added several classes for advanced users.  This can't affect code (if you don't use it)
3678  Made some tiny changes (for N way branching) which should not change anything.
3679  CbcNWay object class - for N way branching this also allows use of CbcConsequence class.
3680  CbcBranchAllDifferent object class - for branching on general integer variables
3681  to stop them having same value so branches are x >= y+1 and x <= y-1.
3682  Added two new Cgl classes - CglAllDifferent which does column fixing (too slowly)
3683  and CglStored which just has a list of cuts which can be activated.
3684  Modified preprocess option to SOS
3685  Version 1.00.02 December 9 2005
3686  Added use of CbcStrategy to do clean preprocessing
3687  Added use of referenceSolver for cleaner repetition of Cbc
3688  Version 1.01.00 February 2 2006
3689  Added first try at Ampl interface
3690*/
Note: See TracBrowser for help on using the repository browser.