source: trunk/Test/CoinSolve.cpp @ 253

Last change on this file since 253 was 253, checked in by forrest, 14 years ago

for verbose help

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