source: trunk/Test/CoinSolve.cpp @ 295

Last change on this file since 295 was 295, checked in by forrest, 15 years ago

for ampl

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