source: trunk/Test/CoinSolve.cpp @ 261

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

for help

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