source: trunk/Test/CoinSolve.cpp @ 254

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

for ampl

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