source: trunk/Test/CoinSolve.cpp @ 252

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

for trysos

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