source: branches/devel/Cbc/src/CoinSolve.cpp @ 465

Last change on this file since 465 was 465, checked in by forrest, 13 years ago

test pseudo shadow prices

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 172.7 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3   
4#include "CbcConfig.h"
5#include "CoinPragma.hpp"
6
7#include <cassert>
8#include <cstdio>
9#include <cmath>
10#include <cfloat>
11#include <cstring>
12#include <iostream>
13
14
15#include "CoinPragma.hpp"
16#include "CoinHelperFunctions.hpp"
17// Same version as CBC
18#define CBCVERSION "1.02.00"
19
20#include "CoinMpsIO.hpp"
21
22#include "ClpFactorization.hpp"
23#include "CoinTime.hpp"
24#include "ClpSimplex.hpp"
25#include "ClpSimplexOther.hpp"
26#include "ClpSolve.hpp"
27#include "ClpPackedMatrix.hpp"
28#include "ClpPlusMinusOneMatrix.hpp"
29#include "ClpNetworkMatrix.hpp"
30#include "ClpDualRowSteepest.hpp"
31#include "ClpDualRowDantzig.hpp"
32#include "ClpLinearObjective.hpp"
33#include "ClpPrimalColumnSteepest.hpp"
34#include "ClpPrimalColumnDantzig.hpp"
35#include "ClpPresolve.hpp"
36#include "CbcOrClpParam.hpp"
37#include "OsiRowCutDebugger.hpp"
38#include "OsiChooseVariable.hpp"
39#ifdef DMALLOC
40#include "dmalloc.h"
41#endif
42#ifdef WSSMP_BARRIER
43#define FOREIGN_BARRIER
44#endif
45#ifdef UFL_BARRIER
46#define FOREIGN_BARRIER
47#endif
48#ifdef TAUCS_BARRIER
49#define FOREIGN_BARRIER
50#endif
51#include "CoinWarmStartBasis.hpp"
52
53#include "OsiSolverInterface.hpp"
54#include "OsiCuts.hpp"
55#include "OsiRowCut.hpp"
56#include "OsiColCut.hpp"
57
58#include "CglPreProcess.hpp"
59#include "CglCutGenerator.hpp"
60#include "CglGomory.hpp"
61#include "CglProbing.hpp"
62#include "CglKnapsackCover.hpp"
63#include "CglRedSplit.hpp"
64#include "CglClique.hpp"
65#include "CglFlowCover.hpp"
66#include "CglMixedIntegerRounding2.hpp"
67#include "CglTwomir.hpp"
68#include "CglDuplicateRow.hpp"
69#include "CglStored.hpp"
70#include "CglLandP.hpp"
71
72#include "CbcModel.hpp"
73#include "CbcHeuristic.hpp"
74#include "CbcHeuristicLocal.hpp"
75#include "CbcHeuristicGreedy.hpp"
76#include "CbcHeuristicFPump.hpp"
77#include "CbcTreeLocal.hpp"
78#include "CbcCompareActual.hpp"
79#include "CbcBranchActual.hpp"
80#include  "CbcOrClpParam.hpp"
81#include  "CbcCutGenerator.hpp"
82#include  "CbcStrategy.hpp"
83
84#include "OsiClpSolverInterface.hpp"
85#ifdef COIN_HAS_ASL
86#include "Cbc_ampl.h"
87static bool usingAmpl=false;
88#endif
89static double totalTime=0.0;
90static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
91static bool maskMatches(const int * starts, char ** masks,
92                        std::string & check);
93static void generateCode(CbcModel * model, const char * fileName,int type,int preProcess);
94#ifdef NDEBUG
95#undef NDEBUG
96#endif
97//#############################################################################
98// To use USERCBC uncomment the following define and add in your fake main program here
99//#define USER_HAS_FAKE_MAIN
100//  Start any fake main program
101#ifdef USER_HAS_FAKE_MAIN
102#endif
103//  End any fake main program
104//#############################################################################
105
106// Allow for interrupts
107// But is this threadsafe ? (so switched off by option)
108
109#include "CoinSignal.hpp"
110static CbcModel * currentBranchModel = NULL;
111
112extern "C" {
113   static void signal_handler(int whichSignal)
114   {
115      if (currentBranchModel!=NULL) 
116         currentBranchModel->setMaximumNodes(0); // stop at next node
117      return;
118   }
119}
120
121int mainTest (int argc, const char *argv[],int algorithm,
122              ClpSimplex empty, bool doPresolve,int switchOff,bool doVector);
123void CbcClpUnitTest (const CbcModel & saveModel);
124int CbcOrClpRead_mode=1;
125FILE * CbcOrClpReadCommand=stdin;
126static bool noPrinting=false;
127static int * analyze(OsiClpSolverInterface * solverMod, int & numberChanged, double & increment,
128                     bool changeInt)
129{
130  OsiSolverInterface * solver = solverMod->clone();
131  if (0) {
132    // just get increment
133    CbcModel model(*solver);
134    model.analyzeObjective();
135    double increment2=model.getCutoffIncrement();
136    printf("initial cutoff increment %g\n",increment2);
137  }
138  const double *objective = solver->getObjCoefficients() ;
139  const double *lower = solver->getColLower() ;
140  const double *upper = solver->getColUpper() ;
141  int numberColumns = solver->getNumCols() ;
142  int numberRows = solver->getNumRows();
143  double direction = solver->getObjSense();
144  int iRow,iColumn;
145
146  // Row copy
147  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
148  const double * elementByRow = matrixByRow.getElements();
149  const int * column = matrixByRow.getIndices();
150  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
151  const int * rowLength = matrixByRow.getVectorLengths();
152
153  // Column copy
154  CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
155  const double * element = matrixByCol.getElements();
156  const int * row = matrixByCol.getIndices();
157  const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
158  const int * columnLength = matrixByCol.getVectorLengths();
159
160  const double * rowLower = solver->getRowLower();
161  const double * rowUpper = solver->getRowUpper();
162
163  char * ignore = new char [numberRows];
164  int * changed = new int[numberColumns];
165  int * which = new int[numberRows];
166  double * changeRhs = new double[numberRows];
167  memset(changeRhs,0,numberRows*sizeof(double));
168  memset(ignore,0,numberRows);
169  numberChanged=0;
170  int numberInteger=0;
171  for (iColumn=0;iColumn<numberColumns;iColumn++) {
172    if (upper[iColumn] > lower[iColumn]+1.0e-8&&solver->isInteger(iColumn)) 
173      numberInteger++;
174  }
175  bool finished=false;
176  while (!finished) {
177    int saveNumberChanged = numberChanged;
178    for (iRow=0;iRow<numberRows;iRow++) {
179      int numberContinuous=0;
180      double value1=0.0,value2=0.0;
181      bool allIntegerCoeff=true;
182      double sumFixed=0.0;
183      int jColumn1=-1,jColumn2=-1;
184      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
185        int jColumn = column[j];
186        double value = elementByRow[j];
187        if (upper[jColumn] > lower[jColumn]+1.0e-8) {
188          if (!solver->isInteger(jColumn)) {
189            if (numberContinuous==0) {
190              jColumn1=jColumn;
191              value1=value;
192            } else {
193              jColumn2=jColumn;
194              value2=value;
195            }
196            numberContinuous++;
197          } else {
198            if (fabs(value-floor(value+0.5))>1.0e-12)
199              allIntegerCoeff=false;
200          }
201        } else {
202          sumFixed += lower[jColumn]*value;
203        }
204      }
205      double low = rowLower[iRow];
206      if (low>-1.0e20) {
207        low -= sumFixed;
208        if (fabs(low-floor(low+0.5))>1.0e-12)
209          allIntegerCoeff=false;
210      }
211      double up = rowUpper[iRow];
212      if (up<1.0e20) {
213        up -= sumFixed;
214        if (fabs(up-floor(up+0.5))>1.0e-12)
215          allIntegerCoeff=false;
216      }
217      if (!allIntegerCoeff)
218        continue; // can't do
219      if (numberContinuous==1) {
220        // see if really integer
221        // This does not allow for complicated cases
222        if (low==up) {
223          if (fabs(value1)>1.0e-3) {
224            value1 = 1.0/value1;
225            if (fabs(value1-floor(value1+0.5))<1.0e-12) {
226              // integer
227              changed[numberChanged++]=jColumn1;
228              solver->setInteger(jColumn1);
229              if (upper[jColumn1]>1.0e20)
230                solver->setColUpper(jColumn1,1.0e20);
231              if (lower[jColumn1]<-1.0e20)
232                solver->setColLower(jColumn1,-1.0e20);
233            }
234          }
235        } else {
236          if (fabs(value1)>1.0e-3) {
237            value1 = 1.0/value1;
238            if (fabs(value1-floor(value1+0.5))<1.0e-12) {
239              // This constraint will not stop it being integer
240              ignore[iRow]=1;
241            }
242          }
243        }
244      } else if (numberContinuous==2) {
245        if (low==up) {
246          /* need general theory - for now just look at 2 cases -
247             1 - +- 1 one in column and just costs i.e. matching objective
248             2 - +- 1 two in column but feeds into G/L row which will try and minimize
249          */
250          if (fabs(value1)==1.0&&value1*value2==-1.0&&!lower[jColumn1]
251              &&!lower[jColumn2]) {
252            int n=0;
253            int i;
254            double objChange=direction*(objective[jColumn1]+objective[jColumn2]);
255            double bound = CoinMin(upper[jColumn1],upper[jColumn2]);
256            bound = CoinMin(bound,1.0e20);
257            for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
258              int jRow = row[i];
259              double value = element[i];
260              if (jRow!=iRow) {
261                which[n++]=jRow;
262                changeRhs[jRow]=value;
263              }
264            }
265            for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
266              int jRow = row[i];
267              double value = element[i];
268              if (jRow!=iRow) {
269                if (!changeRhs[jRow]) {
270                  which[n++]=jRow;
271                  changeRhs[jRow]=value;
272                } else {
273                  changeRhs[jRow]+=value;
274                }
275              }
276            }
277            if (objChange>=0.0) {
278              // see if all rows OK
279              bool good=true;
280              for (i=0;i<n;i++) {
281                int jRow = which[i];
282                double value = changeRhs[jRow];
283                if (value) {
284                  value *= bound;
285                  if (rowLength[jRow]==1) {
286                    if (value>0.0) {
287                      double rhs = rowLower[jRow];
288                      if (rhs>0.0) {
289                        double ratio =rhs/value;
290                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
291                          good=false;
292                      }
293                    } else {
294                      double rhs = rowUpper[jRow];
295                      if (rhs<0.0) {
296                        double ratio =rhs/value;
297                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
298                          good=false;
299                      }
300                    }
301                  } else if (rowLength[jRow]==2) {
302                    if (value>0.0) {
303                      if (rowLower[jRow]>-1.0e20)
304                        good=false;
305                    } else {
306                      if (rowUpper[jRow]<1.0e20)
307                        good=false;
308                    }
309                  } else {
310                    good=false;
311                  }
312                }
313              }
314              if (good) {
315                // both can be integer
316                changed[numberChanged++]=jColumn1;
317                solver->setInteger(jColumn1);
318                if (upper[jColumn1]>1.0e20)
319                  solver->setColUpper(jColumn1,1.0e20);
320                if (lower[jColumn1]<-1.0e20)
321                  solver->setColLower(jColumn1,-1.0e20);
322                changed[numberChanged++]=jColumn2;
323                solver->setInteger(jColumn2);
324                if (upper[jColumn2]>1.0e20)
325                  solver->setColUpper(jColumn2,1.0e20);
326                if (lower[jColumn2]<-1.0e20)
327                  solver->setColLower(jColumn2,-1.0e20);
328              }
329            }
330            // clear
331            for (i=0;i<n;i++) {
332              changeRhs[which[i]]=0.0;
333            }
334          }
335        }
336      }
337    }
338    for (iColumn=0;iColumn<numberColumns;iColumn++) {
339      if (upper[iColumn] > lower[iColumn]+1.0e-8&&!solver->isInteger(iColumn)) {
340        double value;
341        value = upper[iColumn];
342        if (value<1.0e20&&fabs(value-floor(value+0.5))>1.0e-12) 
343          continue;
344        value = lower[iColumn];
345        if (value>-1.0e20&&fabs(value-floor(value+0.5))>1.0e-12) 
346          continue;
347        bool integer=true;
348        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
349          int iRow = row[j];
350          if (!ignore[iRow]) {
351            integer=false;
352            break;
353          }
354        }
355        if (integer) {
356          // integer
357          changed[numberChanged++]=iColumn;
358          solver->setInteger(iColumn);
359          if (upper[iColumn]>1.0e20)
360            solver->setColUpper(iColumn,1.0e20);
361          if (lower[iColumn]<-1.0e20)
362            solver->setColLower(iColumn,-1.0e20);
363        }
364      }
365    }
366    finished = numberChanged==saveNumberChanged;
367  }
368  delete [] which;
369  delete [] changeRhs;
370  delete [] ignore;
371  if (numberInteger&&!noPrinting)
372    printf("%d integer variables",numberInteger);
373  if (changeInt) {
374    if (!noPrinting) {
375      if (numberChanged)
376        printf(" and %d variables made integer\n",numberChanged);
377      else
378        printf("\n");
379    }
380    delete [] ignore;
381    //increment=0.0;
382    if (!numberChanged) {
383      delete [] changed;
384      delete solver;
385      return NULL;
386    } else {
387      for (iColumn=0;iColumn<numberColumns;iColumn++) {
388        if (solver->isInteger(iColumn))
389          solverMod->setInteger(iColumn);
390      }
391      delete solver;
392      return changed;
393    }
394  } else {
395    if (!noPrinting) {
396      if (numberChanged)
397        printf(" and %d variables could be made integer\n",numberChanged);
398      else
399        printf("\n");
400    }
401    // just get increment
402    CbcModel model(*solver);
403    if (noPrinting)
404      model.setLogLevel(0);
405    model.analyzeObjective();
406    double increment2=model.getCutoffIncrement();
407    if (increment2>increment) {
408      if (!noPrinting)
409        printf("cutoff increment increased from %g to %g\n",increment,increment2);
410      increment=increment2;
411    }
412    delete solver;
413    numberChanged=0;
414    delete [] changed;
415    return NULL;
416  }
417}
418static int outDupRow(OsiSolverInterface * solver) 
419{
420  CglDuplicateRow dupCuts(solver);
421  CglTreeInfo info;
422  info.level = 0;
423  info.pass = 0;
424  int numberRows = solver->getNumRows();
425  info.formulation_rows = numberRows;
426  info.inTree = false;
427  info.strengthenRow= NULL;
428  info.pass = 0;
429  OsiCuts cs;
430  dupCuts.generateCuts(*solver,cs,info);
431  const int * duplicate = dupCuts.duplicate();
432  // Get rid of duplicate rows
433  int * which = new int[numberRows]; 
434  int numberDrop=0;
435  for (int iRow=0;iRow<numberRows;iRow++) {
436    if (duplicate[iRow]==-2||duplicate[iRow]>=0) 
437      which[numberDrop++]=iRow;
438  }
439  if (numberDrop) {
440    solver->deleteRows(numberDrop,which);
441  }
442  delete [] which;
443  // see if we have any column cuts
444  int numberColumnCuts = cs.sizeColCuts() ;
445  const double * columnLower = solver->getColLower();
446  const double * columnUpper = solver->getColUpper();
447  for (int k = 0;k<numberColumnCuts;k++) {
448    OsiColCut * thisCut = cs.colCutPtr(k) ;
449    const CoinPackedVector & lbs = thisCut->lbs() ;
450    const CoinPackedVector & ubs = thisCut->ubs() ;
451    int j ;
452    int n ;
453    const int * which ;
454    const double * values ;
455    n = lbs.getNumElements() ;
456    which = lbs.getIndices() ;
457    values = lbs.getElements() ;
458    for (j = 0;j<n;j++) {
459      int iColumn = which[j] ;
460      if (values[j]>columnLower[iColumn]) 
461        solver->setColLower(iColumn,values[j]) ;
462    }
463    n = ubs.getNumElements() ;
464    which = ubs.getIndices() ;
465    values = ubs.getElements() ;
466    for (j = 0;j<n;j++) {
467      int iColumn = which[j] ;
468      if (values[j]<columnUpper[iColumn]) 
469        solver->setColUpper(iColumn,values[j]) ;
470    }
471  }
472  return numberDrop;
473}
474void checkSOS(CbcModel * babModel, const OsiSolverInterface * solver)
475{
476#ifdef COIN_DEVELOP
477  //const double *objective = solver->getObjCoefficients() ;
478  const double *columnLower = solver->getColLower() ;
479  const double * columnUpper = solver->getColUpper() ;
480  const double * solution = solver->getColSolution();
481  //int numberColumns = solver->getNumCols() ;
482  //int numberRows = solver->getNumRows();
483  //double direction = solver->getObjSense();
484  //int iRow,iColumn;
485
486  // Row copy
487  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
488  //const double * elementByRow = matrixByRow.getElements();
489  //const int * column = matrixByRow.getIndices();
490  //const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
491  const int * rowLength = matrixByRow.getVectorLengths();
492
493  // Column copy
494  CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
495  const double * element = matrixByCol.getElements();
496  const int * row = matrixByCol.getIndices();
497  const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
498  const int * columnLength = matrixByCol.getVectorLengths();
499
500  const double * rowLower = solver->getRowLower();
501  const double * rowUpper = solver->getRowUpper();
502  OsiObject ** objects = babModel->objects();
503  int numberObjects = babModel->numberObjects();
504  for (int iObj = 0;iObj<numberObjects;iObj++) {
505    CbcSOS * objSOS =
506      dynamic_cast <CbcSOS *>(objects[iObj]) ;
507    if (objSOS) {
508      int n=objSOS->numberMembers();
509      const int * which = objSOS->members();
510      const double * weight = objSOS->weights();
511      int type = objSOS->sosType();
512      // convexity row?
513      int iColumn;
514      iColumn=which[0];
515      int j;
516      int convex=-1;
517      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
518        int iRow = row[j];
519        double value = element[j];
520        if (rowLower[iRow]==1.0&&rowUpper[iRow]==1.0&&
521            value==1.0) {
522          // possible
523          if (rowLength[iRow]==n) {
524            if (convex==-1)
525              convex=iRow;
526            else
527              convex=-2;
528          }
529        }
530      }
531      printf ("set %d of type %d has %d members - possible convexity row %d\n",
532              iObj,type,n,convex);
533      for (int i=0;i<n;i++) {
534        iColumn = which[i];
535        int convex2=-1;
536        for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
537          int iRow = row[j];
538          if (iRow==convex) {
539            double value = element[j];
540            if (value==1.0) {
541              convex2=iRow;
542            }
543          }
544        }
545        if (convex2<0&&convex>=0) {
546          printf("odd convexity row\n");
547          convex=-2;
548        }
549        printf("col %d has weight %g and value %g, bounds %g %g\n",
550               iColumn,weight[i],solution[iColumn],columnLower[iColumn],
551               columnUpper[iColumn]);
552      }
553    }
554  }
555#endif
556}
557int main (int argc, const char *argv[])
558{
559  /* Note
560     This is meant as a stand-alone executable to do as much of coin as possible.
561     It should only have one solver known to it.
562  */
563  {
564    double time1 = CoinCpuTime(),time2;
565    bool goodModel=false;
566    CoinSighandler_t saveSignal=SIG_DFL;
567    // register signal handler
568    saveSignal = signal(SIGINT,signal_handler);
569    // Set up all non-standard stuff
570    OsiClpSolverInterface solver1;
571    CbcModel model(solver1);
572    CbcModel * babModel = NULL;
573    model.setNumberBeforeTrust(21);
574    int cutPass=-1234567;
575    int tunePreProcess=5;
576    OsiSolverInterface * solver = model.solver();
577    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
578    ClpSimplex * lpSolver = clpSolver->getModelPtr();
579    clpSolver->messageHandler()->setLogLevel(0) ;
580    model.messageHandler()->setLogLevel(1);
581    // For priorities etc
582    int * priorities=NULL;
583    int * branchDirection=NULL;
584    double * pseudoDown=NULL;
585    double * pseudoUp=NULL;
586    double * solutionIn = NULL;
587    int * prioritiesIn = NULL;
588    int numberSOS = 0;
589    int * sosStart = NULL;
590    int * sosIndices = NULL;
591    char * sosType = NULL;
592    double * sosReference = NULL;
593    int * sosPriority=NULL;
594#ifdef COIN_HAS_ASL
595    ampl_info info;
596    memset(&info,0,sizeof(info));
597    if (argc>2&&!strcmp(argv[2],"-AMPL")) {
598      usingAmpl=true;
599      int returnCode = readAmpl(&info,argc,const_cast<char **>(argv));
600      if (returnCode)
601        return returnCode;
602      CbcOrClpRead_mode=2; // so will start with parameters
603      // see if log in list
604      noPrinting=true;
605      for (int i=1;i<info.numberArguments;i++) {
606        if (!strcmp(info.arguments[i],"log")) {
607          if (i<info.numberArguments-1&&atoi(info.arguments[i+1])>0)
608            noPrinting=false;
609          break;
610        }
611      }
612      if (noPrinting) {
613        model.messageHandler()->setLogLevel(0);
614        setCbcOrClpPrinting(false);
615      }
616      if (!noPrinting)
617        printf("%d rows, %d columns and %d elements\n",
618               info.numberRows,info.numberColumns,info.numberElements);
619      solver->loadProblem(info.numberColumns,info.numberRows,info.starts,
620                          info.rows,info.elements,
621                          info.columnLower,info.columnUpper,info.objective,
622                          info.rowLower,info.rowUpper);
623      // If we had a solution use it
624      if (info.primalSolution) {
625        solver->setColSolution(info.primalSolution);
626      }
627      // status
628      if (info.rowStatus) {
629        unsigned char * statusArray = lpSolver->statusArray();
630        int i;
631        for (i=0;i<info.numberColumns;i++)
632          statusArray[i]=(char)info.columnStatus[i];
633        statusArray+=info.numberColumns;
634        for (i=0;i<info.numberRows;i++)
635          statusArray[i]=(char)info.rowStatus[i];
636        CoinWarmStartBasis * basis = lpSolver->getBasis();
637        solver->setWarmStart(basis);
638        delete basis;
639      }
640      freeArrays1(&info);
641      // modify objective if necessary
642      solver->setObjSense(info.direction);
643      solver->setDblParam(OsiObjOffset,info.offset);
644      // Set integer variables
645      for (int i=info.numberColumns-info.numberBinary-info.numberIntegers;
646           i<info.numberColumns;i++)
647        solver->setInteger(i);
648      goodModel=true;
649      // change argc etc
650      argc = info.numberArguments;
651      argv = const_cast<const char **>(info.arguments);
652    }
653#endif   
654    // default action on import
655    int allowImportErrors=0;
656    int keepImportNames=1;
657    int doIdiot=-1;
658    int outputFormat=2;
659    int slpValue=-1;
660    int cppValue=-1;
661    int printOptions=0;
662    int printMode=0;
663    int presolveOptions=0;
664    int substitution=3;
665    int dualize=0;
666    int doCrash=0;
667    int doVector=0;
668    int doSprint=-1;
669    int doScaling=1;
670    // set reasonable defaults
671    int preSolve=5;
672    int preProcess=1;
673    bool useStrategy=false;
674    bool preSolveFile=false;
675   
676    double djFix=1.0e100;
677    double gapRatio=1.0e100;
678    double tightenFactor=0.0;
679    lpSolver->setPerturbation(50);
680    lpSolver->messageHandler()->setPrefix(false);
681    const char dirsep =  CoinFindDirSeparator();
682    std::string directory = (dirsep == '/' ? "./" : ".\\");
683    std::string defaultDirectory = directory;
684    std::string importFile ="";
685    std::string exportFile ="default.mps";
686    std::string importBasisFile ="";
687    std::string importPriorityFile ="";
688    std::string debugFile="";
689    std::string printMask="";
690    double * debugValues = NULL;
691    int numberDebugValues = -1;
692    int basisHasValues=0;
693    std::string exportBasisFile ="default.bas";
694    std::string saveFile ="default.prob";
695    std::string restoreFile ="default.prob";
696    std::string solutionFile ="stdout";
697    std::string solutionSaveFile ="solution.file";
698#define CBCMAXPARAMETERS 200
699    CbcOrClpParam parameters[CBCMAXPARAMETERS];
700    int numberParameters ;
701    establishParams(numberParameters,parameters) ;
702    parameters[whichParam(BASISIN,numberParameters,parameters)].setStringValue(importBasisFile);
703    parameters[whichParam(PRIORITYIN,numberParameters,parameters)].setStringValue(importPriorityFile);
704    parameters[whichParam(BASISOUT,numberParameters,parameters)].setStringValue(exportBasisFile);
705    parameters[whichParam(DEBUG,numberParameters,parameters)].setStringValue(debugFile);
706    parameters[whichParam(PRINTMASK,numberParameters,parameters)].setStringValue(printMask);
707    parameters[whichParam(DIRECTORY,numberParameters,parameters)].setStringValue(directory);
708    parameters[whichParam(DUALBOUND,numberParameters,parameters)].setDoubleValue(lpSolver->dualBound());
709    parameters[whichParam(DUALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->dualTolerance());
710    parameters[whichParam(EXPORT,numberParameters,parameters)].setStringValue(exportFile);
711    parameters[whichParam(IDIOT,numberParameters,parameters)].setIntValue(doIdiot);
712    parameters[whichParam(IMPORT,numberParameters,parameters)].setStringValue(importFile);
713    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].setDoubleValue(1.0e-8);
714    int slog = whichParam(SOLVERLOGLEVEL,numberParameters,parameters);
715    int log = whichParam(LOGLEVEL,numberParameters,parameters);
716    parameters[slog].setIntValue(0);
717    parameters[log].setIntValue(1);
718    parameters[whichParam(MAXFACTOR,numberParameters,parameters)].setIntValue(lpSolver->factorizationFrequency());
719    parameters[whichParam(MAXITERATION,numberParameters,parameters)].setIntValue(lpSolver->maximumIterations());
720    parameters[whichParam(OUTPUTFORMAT,numberParameters,parameters)].setIntValue(outputFormat);
721    parameters[whichParam(PRESOLVEPASS,numberParameters,parameters)].setIntValue(preSolve);
722    parameters[whichParam(PERTVALUE,numberParameters,parameters)].setIntValue(lpSolver->perturbation());
723    parameters[whichParam(PRIMALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->primalTolerance());
724    parameters[whichParam(PRIMALWEIGHT,numberParameters,parameters)].setDoubleValue(lpSolver->infeasibilityCost());
725    parameters[whichParam(RESTORE,numberParameters,parameters)].setStringValue(restoreFile);
726    parameters[whichParam(SAVE,numberParameters,parameters)].setStringValue(saveFile);
727    //parameters[whichParam(TIMELIMIT,numberParameters,parameters)].setDoubleValue(1.0e8);
728    parameters[whichParam(TIMELIMIT_BAB,numberParameters,parameters)].setDoubleValue(1.0e8);
729    parameters[whichParam(SOLUTION,numberParameters,parameters)].setStringValue(solutionFile);
730    parameters[whichParam(SAVESOL,numberParameters,parameters)].setStringValue(solutionSaveFile);
731    parameters[whichParam(SPRINT,numberParameters,parameters)].setIntValue(doSprint);
732    parameters[whichParam(SUBSTITUTION,numberParameters,parameters)].setIntValue(substitution);
733    parameters[whichParam(DUALIZE,numberParameters,parameters)].setIntValue(dualize);
734    model.setNumberBeforeTrust(5);
735    parameters[whichParam(NUMBERBEFORE,numberParameters,parameters)].setIntValue(5);
736    parameters[whichParam(MAXNODES,numberParameters,parameters)].setIntValue(model.getMaximumNodes());
737    model.setNumberStrong(5);
738    parameters[whichParam(STRONGBRANCHING,numberParameters,parameters)].setIntValue(model.numberStrong());
739    parameters[whichParam(INFEASIBILITYWEIGHT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcInfeasibilityWeight));
740    parameters[whichParam(INTEGERTOLERANCE,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcIntegerTolerance));
741    parameters[whichParam(INCREMENT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcCutoffIncrement));
742    // Set up likely cut generators and defaults
743    parameters[whichParam(PREPROCESS,numberParameters,parameters)].setCurrentOption("on");
744    parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].setIntValue(128|64|1);
745    parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].setIntValue(-1);
746    parameters[whichParam(MAXHOTITS,numberParameters,parameters)].setIntValue(100);
747    parameters[whichParam(CUTSSTRATEGY,numberParameters,parameters)].setCurrentOption("on");
748    parameters[whichParam(HEURISTICSTRATEGY,numberParameters,parameters)].setCurrentOption("on");
749    parameters[whichParam(NODESTRATEGY,numberParameters,parameters)].setCurrentOption("fewest");
750    int nodeStrategy=0;
751    int doSOS=1;
752    int verbose=0;
753    CglGomory gomoryGen;
754    // try larger limit
755    gomoryGen.setLimitAtRoot(512);
756    gomoryGen.setLimit(50);
757    // set default action (0=off,1=on,2=root)
758    int gomoryAction=3;
759    parameters[whichParam(GOMORYCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
760
761    CglProbing probingGen;
762    probingGen.setUsingObjective(true);
763    probingGen.setMaxPass(3);
764    probingGen.setMaxPassRoot(3);
765    // Number of unsatisfied variables to look at
766    probingGen.setMaxProbe(10);
767    probingGen.setMaxProbeRoot(50);
768    // How far to follow the consequences
769    probingGen.setMaxLook(10);
770    probingGen.setMaxLookRoot(50);
771    probingGen.setMaxLookRoot(10);
772    // Only look at rows with fewer than this number of elements
773    probingGen.setMaxElements(200);
774    probingGen.setRowCuts(3);
775    // set default action (0=off,1=on,2=root)
776    int probingAction=1;
777    parameters[whichParam(PROBINGCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
778
779    CglKnapsackCover knapsackGen;
780    //knapsackGen.switchOnExpensive();
781    // set default action (0=off,1=on,2=root)
782    int knapsackAction=3;
783    parameters[whichParam(KNAPSACKCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
784
785    CglRedSplit redsplitGen;
786    //redsplitGen.setLimit(100);
787    // set default action (0=off,1=on,2=root)
788    // Off as seems to give some bad cuts
789    int redsplitAction=2;
790    parameters[whichParam(REDSPLITCUTS,numberParameters,parameters)].setCurrentOption("root");
791
792    CglClique cliqueGen(false,true);
793    cliqueGen.setStarCliqueReport(false);
794    cliqueGen.setRowCliqueReport(false);
795    cliqueGen.setMinViolation(0.1);
796    // set default action (0=off,1=on,2=root)
797    int cliqueAction=3;
798    parameters[whichParam(CLIQUECUTS,numberParameters,parameters)].setCurrentOption("ifmove");
799
800    CglMixedIntegerRounding2 mixedGen;
801    // set default action (0=off,1=on,2=root)
802    int mixedAction=3;
803    parameters[whichParam(MIXEDCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
804
805    CglFlowCover flowGen;
806    // set default action (0=off,1=on,2=root)
807    int flowAction=3;
808    parameters[whichParam(FLOWCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
809
810    CglTwomir twomirGen;
811    twomirGen.setMaxElements(250);
812    // set default action (0=off,1=on,2=root)
813    int twomirAction=2;
814    parameters[whichParam(TWOMIRCUTS,numberParameters,parameters)].setCurrentOption("root");
815    CglLandP landpGen;
816    // set default action (0=off,1=on,2=root)
817    int landpAction=0;
818    parameters[whichParam(LANDPCUTS,numberParameters,parameters)].setCurrentOption("off");
819    // Stored cuts
820    bool storedCuts = false;
821
822    bool useRounding=true;
823    parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption("on");
824    bool useFpump=true;
825    parameters[whichParam(FPUMP,numberParameters,parameters)].setCurrentOption("on");
826    bool useGreedy=true;
827    parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption("on");
828    bool useCombine=true;
829    parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption("on");
830    bool useLocalTree=false;
831    parameters[whichParam(COSTSTRATEGY,numberParameters,parameters)].setCurrentOption("off");
832    int useCosts=0;
833    // don't use input solution
834    int useSolution=0;
835   
836    // total number of commands read
837    int numberGoodCommands=0;
838    // Set false if user does anything advanced
839    bool defaultSettings=true;
840
841    // Hidden stuff for barrier
842    int choleskyType = 0;
843    int gamma=0;
844    int scaleBarrier=0;
845    int doKKT=0;
846    int crossover=2; // do crossover unless quadratic
847    // For names
848    int lengthName = 0;
849    std::vector<std::string> rowNames;
850    std::vector<std::string> columnNames;
851   
852    std::string field;
853    if (!noPrinting) {
854      std::cout<<"Coin Cbc and Clp Solver version "<<CBCVERSION
855               <<", build "<<__DATE__<<std::endl;
856      // Print command line
857      if (argc>1) {
858        printf("command line - ");
859        for (int i=0;i<argc;i++)
860          printf("%s ",argv[i]);
861        printf("\n");
862      }
863    }
864    while (1) {
865      // next command
866      field=CoinReadGetCommand(argc,argv);
867      // exit if null or similar
868      if (!field.length()) {
869        if (numberGoodCommands==1&&goodModel) {
870          // we just had file name - do branch and bound
871          field="branch";
872        } else if (!numberGoodCommands) {
873          // let's give the sucker a hint
874          std::cout
875            <<"CoinSolver takes input from arguments ( - switches to stdin)"
876            <<std::endl
877            <<"Enter ? for list of commands or help"<<std::endl;
878          field="-";
879        } else {
880          break;
881        }
882      }
883     
884      // see if ? at end
885      int numberQuery=0;
886      if (field!="?"&&field!="???") {
887        int length = field.length();
888        int i;
889        for (i=length-1;i>0;i--) {
890          if (field[i]=='?') 
891            numberQuery++;
892          else
893            break;
894        }
895        field=field.substr(0,length-numberQuery);
896      }
897      // find out if valid command
898      int iParam;
899      int numberMatches=0;
900      int firstMatch=-1;
901      for ( iParam=0; iParam<numberParameters; iParam++ ) {
902        int match = parameters[iParam].matches(field);
903        if (match==1) {
904          numberMatches = 1;
905          firstMatch=iParam;
906          break;
907        } else {
908          if (match&&firstMatch<0)
909            firstMatch=iParam;
910          numberMatches += match>>1;
911        }
912      }
913      if (iParam<numberParameters&&!numberQuery) {
914        // found
915        CbcOrClpParam found = parameters[iParam];
916        CbcOrClpParameterType type = found.type();
917        int valid;
918        numberGoodCommands++;
919        if (type==BAB&&goodModel) {
920          // check if any integers
921#ifdef COIN_HAS_ASL
922          if (info.numberSos&&doSOS&&usingAmpl) {
923            // SOS
924            numberSOS = info.numberSos;
925          }
926#endif
927          if (!lpSolver->integerInformation()&&!numberSOS&&
928              !clpSolver->numberSOS())
929            type=DUALSIMPLEX;
930        }
931        if (type==GENERALQUERY) {
932          bool evenHidden=false;
933          if ((verbose&8)!=0) {
934            // even hidden
935            evenHidden = true;
936            verbose &= ~8;
937          }
938#ifdef COIN_HAS_ASL
939          if (verbose<4&&usingAmpl)
940            verbose +=4;
941#endif
942          if (verbose<4) {
943            std::cout<<"In argument list keywords have leading - "
944              ", -stdin or just - switches to stdin"<<std::endl;
945            std::cout<<"One command per line (and no -)"<<std::endl;
946            std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
947            std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
948            std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
949            std::cout<<"abcd value sets value"<<std::endl;
950            std::cout<<"Commands are:"<<std::endl;
951          } else {
952            std::cout<<"Cbc options are set within AMPL with commands like:"<<std::endl<<std::endl;
953            std::cout<<"         option cbc_options \"cuts=root log=2 feas=on slog=1\""<<std::endl<<std::endl;
954            std::cout<<"only maximize, dual, primal, help and quit are recognized without ="<<std::endl;
955          }
956          int maxAcross=5;
957          if ((verbose%4)!=0)
958            maxAcross=1;
959          int limits[]={1,51,101,151,201,251,301,351,401};
960          std::vector<std::string> types;
961          types.push_back("Double parameters:");
962          types.push_back("Branch and Cut double parameters:");
963          types.push_back("Integer parameters:");
964          types.push_back("Branch and Cut integer parameters:");
965          types.push_back("Keyword parameters:");
966          types.push_back("Branch and Cut keyword parameters:");
967          types.push_back("Actions or string parameters:");
968          types.push_back("Branch and Cut actions:");
969          int iType;
970          for (iType=0;iType<8;iType++) {
971            int across=0;
972            if ((verbose%4)!=0)
973              std::cout<<std::endl;
974            std::cout<<types[iType]<<std::endl;
975            if ((verbose&2)!=0)
976              std::cout<<std::endl;
977            for ( iParam=0; iParam<numberParameters; iParam++ ) {
978              int type = parameters[iParam].type();
979              if ((parameters[iParam].displayThis()||evenHidden)&&
980                  type>=limits[iType]
981                  &&type<limits[iType+1]) {
982                // but skip if not useful for ampl (and in ampl mode)
983                if (verbose>=4&&(parameters[iParam].whereUsed()&4)==0)
984                  continue;
985                if (!across) {
986                  if ((verbose&2)==0) 
987                    std::cout<<"  ";
988                  else
989                    std::cout<<"Command ";
990                }
991                std::cout<<parameters[iParam].matchName()<<"  ";
992                across++;
993                if (across==maxAcross) {
994                  across=0;
995                  if ((verbose%4)!=0) {
996                    // put out description as well
997                    if ((verbose&1)!=0) 
998                      std::cout<<parameters[iParam].shortHelp();
999                    std::cout<<std::endl;
1000                    if ((verbose&2)!=0) {
1001                      std::cout<<"---- description"<<std::endl;
1002                      parameters[iParam].printLongHelp();
1003                      std::cout<<"----"<<std::endl<<std::endl;
1004                    }
1005                  } else {
1006                    std::cout<<std::endl;
1007                  }
1008                }
1009              }
1010            }
1011            if (across)
1012              std::cout<<std::endl;
1013          }
1014        } else if (type==FULLGENERALQUERY) {
1015          std::cout<<"Full list of commands is:"<<std::endl;
1016          int maxAcross=5;
1017          int limits[]={1,51,101,151,201,251,301,351,401};
1018          std::vector<std::string> types;
1019          types.push_back("Double parameters:");
1020          types.push_back("Branch and Cut double parameters:");
1021          types.push_back("Integer parameters:");
1022          types.push_back("Branch and Cut integer parameters:");
1023          types.push_back("Keyword parameters:");
1024          types.push_back("Branch and Cut keyword parameters:");
1025          types.push_back("Actions or string parameters:");
1026          types.push_back("Branch and Cut actions:");
1027          int iType;
1028          for (iType=0;iType<8;iType++) {
1029            int across=0;
1030            std::cout<<types[iType]<<"  ";
1031            for ( iParam=0; iParam<numberParameters; iParam++ ) {
1032              int type = parameters[iParam].type();
1033              if (type>=limits[iType]
1034                  &&type<limits[iType+1]) {
1035                if (!across)
1036                  std::cout<<"  ";
1037                std::cout<<parameters[iParam].matchName()<<"  ";
1038                across++;
1039                if (across==maxAcross) {
1040                  std::cout<<std::endl;
1041                  across=0;
1042                }
1043              }
1044            }
1045            if (across)
1046              std::cout<<std::endl;
1047          }
1048        } else if (type<101) {
1049          // get next field as double
1050          double value = CoinReadGetDoubleField(argc,argv,&valid);
1051          if (!valid) {
1052            if (type<51) {
1053              parameters[iParam].setDoubleParameter(lpSolver,value);
1054            } else if (type<81) {
1055              parameters[iParam].setDoubleParameter(model,value);
1056            } else {
1057              parameters[iParam].setDoubleParameter(lpSolver,value);
1058              switch(type) {
1059              case DJFIX:
1060                djFix=value;
1061                preSolve=5;
1062                defaultSettings=false; // user knows what she is doing
1063                break;
1064              case GAPRATIO:
1065                gapRatio=value;
1066                break;
1067              case TIGHTENFACTOR:
1068                tightenFactor=value;
1069                defaultSettings=false; // user knows what she is doing
1070                break;
1071              default:
1072                abort();
1073              }
1074            }
1075          } else if (valid==1) {
1076            abort();
1077          } else {
1078            std::cout<<parameters[iParam].name()<<" has value "<<
1079              parameters[iParam].doubleValue()<<std::endl;
1080          }
1081        } else if (type<201) {
1082          // get next field as int
1083          int value = CoinReadGetIntField(argc,argv,&valid);
1084          if (!valid) {
1085            if (type<151) {
1086              if (parameters[iParam].type()==PRESOLVEPASS)
1087                preSolve = value;
1088              else if (parameters[iParam].type()==IDIOT)
1089                doIdiot = value;
1090              else if (parameters[iParam].type()==SPRINT)
1091                doSprint = value;
1092              else if (parameters[iParam].type()==OUTPUTFORMAT)
1093                outputFormat = value;
1094              else if (parameters[iParam].type()==SLPVALUE)
1095                slpValue = value;
1096              else if (parameters[iParam].type()==CPP)
1097                cppValue = value;
1098              else if (parameters[iParam].type()==PRESOLVEOPTIONS)
1099                presolveOptions = value;
1100              else if (parameters[iParam].type()==PRINTOPTIONS)
1101                printOptions = value;
1102              else if (parameters[iParam].type()==SUBSTITUTION)
1103                substitution = value;
1104              else if (parameters[iParam].type()==DUALIZE)
1105                dualize = value;
1106              else if (parameters[iParam].type()==CUTPASS)
1107                cutPass = value;
1108              else if (parameters[iParam].type()==PROCESSTUNE)
1109                tunePreProcess = value;
1110              else if (parameters[iParam].type()==VERBOSE)
1111                verbose = value;
1112              else if (parameters[iParam].type()==FPUMPITS)
1113                { useFpump = true;parameters[iParam].setIntValue(value);}
1114              parameters[iParam].setIntParameter(lpSolver,value);
1115            } else {
1116              parameters[iParam].setIntParameter(model,value);
1117            }
1118          } else if (valid==1) {
1119            abort();
1120          } else {
1121            std::cout<<parameters[iParam].name()<<" has value "<<
1122              parameters[iParam].intValue()<<std::endl;
1123          }
1124        } else if (type<301) {
1125          // one of several strings
1126          std::string value = CoinReadGetString(argc,argv);
1127          int action = parameters[iParam].parameterOption(value);
1128          if (action<0) {
1129            if (value!="EOL") {
1130              // no match
1131              parameters[iParam].printOptions();
1132            } else {
1133              // print current value
1134              std::cout<<parameters[iParam].name()<<" has value "<<
1135                parameters[iParam].currentOption()<<std::endl;
1136            }
1137          } else {
1138            parameters[iParam].setCurrentOption(action,!noPrinting);
1139            // for now hard wired
1140            switch (type) {
1141            case DIRECTION:
1142              if (action==0)
1143                lpSolver->setOptimizationDirection(1);
1144              else if (action==1)
1145                lpSolver->setOptimizationDirection(-1);
1146              else
1147                lpSolver->setOptimizationDirection(0);
1148              break;
1149            case DUALPIVOT:
1150              if (action==0) {
1151                ClpDualRowSteepest steep(3);
1152                lpSolver->setDualRowPivotAlgorithm(steep);
1153              } else if (action==1) {
1154                //ClpDualRowDantzig dantzig;
1155                ClpDualRowSteepest dantzig(5);
1156                lpSolver->setDualRowPivotAlgorithm(dantzig);
1157              } else if (action==2) {
1158                // partial steep
1159                ClpDualRowSteepest steep(2);
1160                lpSolver->setDualRowPivotAlgorithm(steep);
1161              } else {
1162                ClpDualRowSteepest steep;
1163                lpSolver->setDualRowPivotAlgorithm(steep);
1164              }
1165              break;
1166            case PRIMALPIVOT:
1167              if (action==0) {
1168                ClpPrimalColumnSteepest steep(3);
1169                lpSolver->setPrimalColumnPivotAlgorithm(steep);
1170              } else if (action==1) {
1171                ClpPrimalColumnSteepest steep(0);
1172                lpSolver->setPrimalColumnPivotAlgorithm(steep);
1173              } else if (action==2) {
1174                ClpPrimalColumnDantzig dantzig;
1175                lpSolver->setPrimalColumnPivotAlgorithm(dantzig);
1176              } else if (action==3) {
1177                ClpPrimalColumnSteepest steep(2);
1178                lpSolver->setPrimalColumnPivotAlgorithm(steep);
1179              } else if (action==4) {
1180                ClpPrimalColumnSteepest steep(1);
1181                lpSolver->setPrimalColumnPivotAlgorithm(steep);
1182              } else if (action==5) {
1183                ClpPrimalColumnSteepest steep(4);
1184                lpSolver->setPrimalColumnPivotAlgorithm(steep);
1185              } else if (action==6) {
1186                ClpPrimalColumnSteepest steep(10);
1187                lpSolver->setPrimalColumnPivotAlgorithm(steep);
1188              }
1189              break;
1190            case SCALING:
1191              lpSolver->scaling(action);
1192              solver->setHintParam(OsiDoScale,action!=0,OsiHintTry);
1193              doScaling = 1-action;
1194              break;
1195            case AUTOSCALE:
1196              lpSolver->setAutomaticScaling(action!=0);
1197              break;
1198            case SPARSEFACTOR:
1199              lpSolver->setSparseFactorization((1-action)!=0);
1200              break;
1201            case BIASLU:
1202              lpSolver->factorization()->setBiasLU(action);
1203              break;
1204            case PERTURBATION:
1205              if (action==0)
1206                lpSolver->setPerturbation(50);
1207              else
1208                lpSolver->setPerturbation(100);
1209              break;
1210            case ERRORSALLOWED:
1211              allowImportErrors = action;
1212              break;
1213            case INTPRINT:
1214              printMode=action;
1215              break;
1216              //case ALGORITHM:
1217              //algorithm  = action;
1218              //defaultSettings=false; // user knows what she is doing
1219              //abort();
1220              //break;
1221            case KEEPNAMES:
1222              keepImportNames = 1-action;
1223              break;
1224            case PRESOLVE:
1225              if (action==0)
1226                preSolve = 5;
1227              else if (action==1)
1228                preSolve=0;
1229              else if (action==2)
1230                preSolve=10;
1231              else
1232                preSolveFile=true;
1233              break;
1234            case PFI:
1235              lpSolver->factorization()->setForrestTomlin(action==0);
1236              break;
1237            case CRASH:
1238              doCrash=action;
1239              break;
1240            case VECTOR:
1241              doVector=action;
1242              break;
1243            case MESSAGES:
1244              lpSolver->messageHandler()->setPrefix(action!=0);
1245              break;
1246            case CHOLESKY:
1247              choleskyType = action;
1248              break;
1249            case GAMMA:
1250              gamma=action;
1251              break;
1252            case BARRIERSCALE:
1253              scaleBarrier=action;
1254              break;
1255            case KKT:
1256              doKKT=action;
1257              break;
1258            case CROSSOVER:
1259              crossover=action;
1260              break;
1261            case SOS:
1262              doSOS=action;
1263              break;
1264            case GOMORYCUTS:
1265              defaultSettings=false; // user knows what she is doing
1266              gomoryAction = action;
1267              break;
1268            case PROBINGCUTS:
1269              defaultSettings=false; // user knows what she is doing
1270              probingAction = action;
1271              break;
1272            case KNAPSACKCUTS:
1273              defaultSettings=false; // user knows what she is doing
1274              knapsackAction = action;
1275              break;
1276            case REDSPLITCUTS:
1277              defaultSettings=false; // user knows what she is doing
1278              redsplitAction = action;
1279              break;
1280            case CLIQUECUTS:
1281              defaultSettings=false; // user knows what she is doing
1282              cliqueAction = action;
1283              break;
1284            case FLOWCUTS:
1285              defaultSettings=false; // user knows what she is doing
1286              flowAction = action;
1287              break;
1288            case MIXEDCUTS:
1289              defaultSettings=false; // user knows what she is doing
1290              mixedAction = action;
1291              break;
1292            case TWOMIRCUTS:
1293              defaultSettings=false; // user knows what she is doing
1294              twomirAction = action;
1295              break;
1296            case LANDPCUTS:
1297              defaultSettings=false; // user knows what she is doing
1298              landpAction = action;
1299              break;
1300            case ROUNDING:
1301              defaultSettings=false; // user knows what she is doing
1302              useRounding = action;
1303              break;
1304            case FPUMP:
1305              defaultSettings=false; // user knows what she is doing
1306              useFpump=action;
1307              break;
1308            case CUTSSTRATEGY:
1309              gomoryAction = action;
1310              probingAction = action;
1311              knapsackAction = action;
1312              cliqueAction = action;
1313              flowAction = action;
1314              mixedAction = action;
1315              twomirAction = action;
1316              //landpAction = action;
1317              parameters[whichParam(GOMORYCUTS,numberParameters,parameters)].setCurrentOption(action);
1318              parameters[whichParam(PROBINGCUTS,numberParameters,parameters)].setCurrentOption(action);
1319              parameters[whichParam(KNAPSACKCUTS,numberParameters,parameters)].setCurrentOption(action);
1320              if (!action) {
1321                redsplitAction = action;
1322                parameters[whichParam(REDSPLITCUTS,numberParameters,parameters)].setCurrentOption(action);
1323              }
1324              parameters[whichParam(CLIQUECUTS,numberParameters,parameters)].setCurrentOption(action);
1325              parameters[whichParam(FLOWCUTS,numberParameters,parameters)].setCurrentOption(action);
1326              parameters[whichParam(MIXEDCUTS,numberParameters,parameters)].setCurrentOption(action);
1327              parameters[whichParam(TWOMIRCUTS,numberParameters,parameters)].setCurrentOption(action);
1328              if (!action) {
1329                landpAction = action;
1330                parameters[whichParam(LANDPCUTS,numberParameters,parameters)].setCurrentOption(action);
1331              }
1332              break;
1333            case HEURISTICSTRATEGY:
1334              useRounding = action;
1335              useGreedy = action;
1336              useCombine = action;
1337              //useLocalTree = action;
1338              useFpump=action;
1339              parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption(action);
1340              parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption(action);
1341              parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption(action);
1342              //parameters[whichParam(LOCALTREE,numberParameters,parameters)].setCurrentOption(action);
1343              parameters[whichParam(FPUMP,numberParameters,parameters)].setCurrentOption(action);
1344              break;
1345            case GREEDY:
1346              defaultSettings=false; // user knows what she is doing
1347              useGreedy = action;
1348              break;
1349            case COMBINE:
1350              defaultSettings=false; // user knows what she is doing
1351              useCombine = action;
1352              break;
1353            case LOCALTREE:
1354              defaultSettings=false; // user knows what she is doing
1355              useLocalTree = action;
1356              break;
1357            case COSTSTRATEGY:
1358              useCosts=action;
1359              break;
1360            case NODESTRATEGY:
1361              nodeStrategy=action;
1362              break;
1363            case PREPROCESS:
1364              preProcess = action;
1365              break;
1366            case USESOLUTION:
1367              useSolution = action;
1368              break;
1369            default:
1370              abort();
1371            }
1372          }
1373        } else {
1374          // action
1375          if (type==EXIT) {
1376#ifdef COIN_HAS_ASL
1377            if(usingAmpl) {
1378              if (info.numberIntegers||info.numberBinary) {
1379                // integer
1380              } else {
1381                // linear
1382              }
1383              writeAmpl(&info);
1384              freeArrays2(&info);
1385              freeArgs(&info);
1386            }
1387#endif
1388            break; // stop all
1389          }
1390          switch (type) {
1391          case DUALSIMPLEX:
1392          case PRIMALSIMPLEX:
1393          case SOLVECONTINUOUS:
1394          case BARRIER:
1395            if (goodModel) {
1396              double objScale = 
1397                parameters[whichParam(OBJSCALE2,numberParameters,parameters)].doubleValue();
1398              if (objScale!=1.0) {
1399                int iColumn;
1400                int numberColumns=lpSolver->numberColumns();
1401                double * dualColumnSolution = 
1402                  lpSolver->dualColumnSolution();
1403                ClpObjective * obj = lpSolver->objectiveAsObject();
1404                assert(dynamic_cast<ClpLinearObjective *> (obj));
1405                double offset;
1406                double * objective = obj->gradient(NULL,NULL,offset,true);
1407                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1408                  dualColumnSolution[iColumn] *= objScale;
1409                  objective[iColumn] *= objScale;;
1410                }
1411                int iRow;
1412                int numberRows=lpSolver->numberRows();
1413                double * dualRowSolution = 
1414                  lpSolver->dualRowSolution();
1415                for (iRow=0;iRow<numberRows;iRow++) 
1416                  dualRowSolution[iRow] *= objScale;
1417                lpSolver->setObjectiveOffset(objScale*lpSolver->objectiveOffset());
1418              }
1419              ClpSolve::SolveType method;
1420              ClpSolve::PresolveType presolveType;
1421              ClpSimplex * model2 = lpSolver;
1422              if (dualize) {
1423                model2 = ((ClpSimplexOther *) model2)->dualOfModel();
1424                printf("Dual of model has %d rows and %d columns\n",
1425                       model2->numberRows(),model2->numberColumns());
1426                model2->setOptimizationDirection(1.0);
1427              }
1428              if (noPrinting)
1429                lpSolver->setLogLevel(0);
1430              ClpSolve solveOptions;
1431              solveOptions.setPresolveActions(presolveOptions);
1432              solveOptions.setSubstitution(substitution);
1433              if (preSolve!=5&&preSolve) {
1434                presolveType=ClpSolve::presolveNumber;
1435                if (preSolve<0) {
1436                  preSolve = - preSolve;
1437                  if (preSolve<=100) {
1438                    presolveType=ClpSolve::presolveNumber;
1439                    printf("Doing %d presolve passes - picking up non-costed slacks\n",
1440                           preSolve);
1441                    solveOptions.setDoSingletonColumn(true);
1442                  } else {
1443                    preSolve -=100;
1444                    presolveType=ClpSolve::presolveNumberCost;
1445                    printf("Doing %d presolve passes - picking up costed slacks\n",
1446                           preSolve);
1447                  }
1448                } 
1449              } else if (preSolve) {
1450                presolveType=ClpSolve::presolveOn;
1451              } else {
1452                presolveType=ClpSolve::presolveOff;
1453              }
1454              solveOptions.setPresolveType(presolveType,preSolve);
1455              if (type==DUALSIMPLEX||type==SOLVECONTINUOUS) {
1456                method=ClpSolve::useDual;
1457              } else if (type==PRIMALSIMPLEX) {
1458                method=ClpSolve::usePrimalorSprint;
1459              } else {
1460                method = ClpSolve::useBarrier;
1461                if (crossover==1) {
1462                  method=ClpSolve::useBarrierNoCross;
1463                } else if (crossover==2) {
1464                  ClpObjective * obj = lpSolver->objectiveAsObject();
1465                  if (obj->type()>1) {
1466                    method=ClpSolve::useBarrierNoCross;
1467                    presolveType=ClpSolve::presolveOff;
1468                    solveOptions.setPresolveType(presolveType,preSolve);
1469                  } 
1470                }
1471              }
1472              solveOptions.setSolveType(method);
1473              if(preSolveFile)
1474                presolveOptions |= 0x40000000;
1475              solveOptions.setSpecialOption(4,presolveOptions);
1476              solveOptions.setSpecialOption(5,printOptions);
1477              if (doVector) {
1478                ClpMatrixBase * matrix = lpSolver->clpMatrix();
1479                if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
1480                  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1481                  clpMatrix->makeSpecialColumnCopy();
1482                }
1483              }
1484              if (method==ClpSolve::useDual) {
1485                // dual
1486                if (doCrash)
1487                  solveOptions.setSpecialOption(0,1,doCrash); // crash
1488                else if (doIdiot)
1489                  solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
1490              } else if (method==ClpSolve::usePrimalorSprint) {
1491                // primal
1492                // if slp turn everything off
1493                if (slpValue>0) {
1494                  doCrash=false;
1495                  doSprint=0;
1496                  doIdiot=-1;
1497                  solveOptions.setSpecialOption(1,10,slpValue); // slp
1498                  method=ClpSolve::usePrimal;
1499                }
1500                if (doCrash) {
1501                  solveOptions.setSpecialOption(1,1,doCrash); // crash
1502                } else if (doSprint>0) {
1503                  // sprint overrides idiot
1504                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
1505                } else if (doIdiot>0) {
1506                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
1507                } else if (slpValue<=0) {
1508                  if (doIdiot==0) {
1509                    if (doSprint==0)
1510                      solveOptions.setSpecialOption(1,4); // all slack
1511                    else
1512                      solveOptions.setSpecialOption(1,9); // all slack or sprint
1513                  } else {
1514                    if (doSprint==0)
1515                      solveOptions.setSpecialOption(1,8); // all slack or idiot
1516                    else
1517                      solveOptions.setSpecialOption(1,7); // initiative
1518                  }
1519                }
1520                if (basisHasValues==-1)
1521                  solveOptions.setSpecialOption(1,11); // switch off values
1522              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
1523                int barrierOptions = choleskyType;
1524                if (scaleBarrier)
1525                  barrierOptions |= 8;
1526                if (doKKT)
1527                  barrierOptions |= 16;
1528                if (gamma)
1529                  barrierOptions |= 32*gamma;
1530                if (crossover==3) 
1531                  barrierOptions |= 256; // try presolve in crossover
1532                solveOptions.setSpecialOption(4,barrierOptions);
1533              }
1534              model2->initialSolve(solveOptions);
1535              basisHasValues=1;
1536              if (dualize) {
1537                int returnCode=((ClpSimplexOther *) lpSolver)->restoreFromDual(model2);
1538                delete model2;
1539                if (returnCode)
1540                  lpSolver->primal(1);
1541                model2=lpSolver;
1542              }
1543#ifdef COIN_HAS_ASL
1544              if (usingAmpl) {
1545                double value = model2->getObjValue()*model2->getObjSense();
1546                char buf[300];
1547                int pos=0;
1548                int iStat = model2->status();
1549                if (iStat==0) {
1550                  pos += sprintf(buf+pos,"optimal," );
1551                } else if (iStat==1) {
1552                  // infeasible
1553                  pos += sprintf(buf+pos,"infeasible,");
1554                } else if (iStat==2) {
1555                  // unbounded
1556                  pos += sprintf(buf+pos,"unbounded,");
1557                } else if (iStat==3) {
1558                  pos += sprintf(buf+pos,"stopped on iterations or time,");
1559                } else if (iStat==4) {
1560                  iStat = 7;
1561                  pos += sprintf(buf+pos,"stopped on difficulties,");
1562                } else if (iStat==5) {
1563                  iStat = 3;
1564                  pos += sprintf(buf+pos,"stopped on ctrl-c,");
1565                } else {
1566                  pos += sprintf(buf+pos,"status unknown,");
1567                  iStat=6;
1568                }
1569                info.problemStatus=iStat;
1570                info.objValue = value;
1571                pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
1572                               value);
1573                sprintf(buf+pos,"\n%d iterations",
1574                        model2->getIterationCount());
1575                free(info.primalSolution);
1576                int numberColumns=model2->numberColumns();
1577                info.primalSolution = (double *) malloc(numberColumns*sizeof(double));
1578                CoinCopyN(model2->primalColumnSolution(),numberColumns,info.primalSolution);
1579                int numberRows = model2->numberRows();
1580                free(info.dualSolution);
1581                info.dualSolution = (double *) malloc(numberRows*sizeof(double));
1582                CoinCopyN(model2->dualRowSolution(),numberRows,info.dualSolution);
1583                CoinWarmStartBasis * basis = model2->getBasis();
1584                free(info.rowStatus);
1585                info.rowStatus = (int *) malloc(numberRows*sizeof(int));
1586                free(info.columnStatus);
1587                info.columnStatus = (int *) malloc(numberColumns*sizeof(int));
1588                // Put basis in
1589                int i;
1590                // free,basic,ub,lb are 0,1,2,3
1591                for (i=0;i<numberRows;i++) {
1592                  CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
1593                  info.rowStatus[i]=status;
1594                }
1595                for (i=0;i<numberColumns;i++) {
1596                  CoinWarmStartBasis::Status status = basis->getStructStatus(i);
1597                  info.columnStatus[i]=status;
1598                }
1599                // put buffer into info
1600                strcpy(info.buffer,buf);
1601                delete basis;
1602              }
1603#endif
1604            } else {
1605              std::cout<<"** Current model not valid"<<std::endl;
1606            }
1607            break;
1608          case STATISTICS:
1609            if (goodModel) {
1610              // If presolve on look at presolved
1611              bool deleteModel2=false;
1612              ClpSimplex * model2 = lpSolver;
1613              if (preSolve) {
1614                ClpPresolve pinfo;
1615                int presolveOptions2 = presolveOptions&~0x40000000;
1616                if ((presolveOptions2&0xffff)!=0)
1617                  pinfo.setPresolveActions(presolveOptions2);
1618                pinfo.setSubstitution(substitution);
1619                if ((printOptions&1)!=0)
1620                  pinfo.statistics();
1621                double presolveTolerance = 
1622                  parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1623                model2 = 
1624                  pinfo.presolvedModel(*lpSolver,presolveTolerance,
1625                                       true,preSolve);
1626                if (model2) {
1627                  printf("Statistics for presolved model\n");
1628                  deleteModel2=true;
1629                } else {
1630                  printf("Presolved model looks infeasible - will use unpresolved\n");
1631                  model2 = lpSolver;
1632                }
1633              } else {
1634                printf("Statistics for unpresolved model\n");
1635                model2 =  lpSolver;
1636              }
1637              statistics(lpSolver,model2);
1638              if (deleteModel2)
1639                delete model2;
1640            } else {
1641              std::cout<<"** Current model not valid"<<std::endl;
1642            }
1643            break;
1644          case TIGHTEN:
1645            if (goodModel) {
1646              int numberInfeasibilities = lpSolver->tightenPrimalBounds();
1647              if (numberInfeasibilities)
1648                std::cout<<"** Analysis indicates model infeasible"<<std::endl;
1649            } else {
1650              std::cout<<"** Current model not valid"<<std::endl;
1651            }
1652            break;
1653          case PLUSMINUS:
1654            if (goodModel) {
1655              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
1656              ClpPackedMatrix* clpMatrix =
1657                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1658              if (clpMatrix) {
1659                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
1660                if (newMatrix->getIndices()) {
1661                  lpSolver->replaceMatrix(newMatrix);
1662                  delete saveMatrix;
1663                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
1664                } else {
1665                  std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
1666                }
1667              } else {
1668                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
1669              }
1670            } else {
1671              std::cout<<"** Current model not valid"<<std::endl;
1672            }
1673            break;
1674          case OUTDUPROWS:
1675            if (goodModel) {
1676              int numberRows = clpSolver->getNumRows();
1677              //int nOut = outDupRow(clpSolver);
1678              CglDuplicateRow dupcuts(clpSolver);
1679              storedCuts = dupcuts.outDuplicates(clpSolver)!=0;
1680              int nOut = numberRows-clpSolver->getNumRows();
1681              if (nOut&&!noPrinting)
1682                printf("%d rows eliminated\n",nOut);
1683            } else {
1684              std::cout<<"** Current model not valid"<<std::endl;
1685            }
1686            break;
1687          case NETWORK:
1688            if (goodModel) {
1689              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
1690              ClpPackedMatrix* clpMatrix =
1691                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1692              if (clpMatrix) {
1693                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
1694                if (newMatrix->getIndices()) {
1695                  lpSolver->replaceMatrix(newMatrix);
1696                  delete saveMatrix;
1697                  std::cout<<"Matrix converted to network matrix"<<std::endl;
1698                } else {
1699                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
1700                }
1701              } else {
1702                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
1703              }
1704            } else {
1705              std::cout<<"** Current model not valid"<<std::endl;
1706            }
1707            break;
1708          case MIPLIB:
1709            // User can set options - main difference is lack of model and CglPreProcess
1710            goodModel=true;
1711/*
1712  Run branch-and-cut. First set a few options -- node comparison, scaling. If
1713  the solver is Clp, consider running some presolve code (not yet converted
1714  this to generic OSI) with branch-and-cut. If presolve is disabled, or the
1715  solver is not Clp, simply run branch-and-cut. Print elapsed time at the end.
1716*/
1717          case BAB: // branchAndBound
1718          case STRENGTHEN:
1719            if (goodModel) {
1720              bool miplib = type==MIPLIB;
1721              int logLevel = parameters[slog].intValue();
1722              // Reduce printout
1723              if (logLevel<=1)
1724                model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
1725              // Don't switch off all output
1726              {
1727                OsiSolverInterface * solver = model.solver();
1728                OsiClpSolverInterface * si =
1729                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
1730                assert (si != NULL);
1731                si->setSpecialOptions(0x40000000);
1732              }
1733              if (!miplib) {
1734                model.initialSolve();
1735                OsiSolverInterface * solver = model.solver();
1736                OsiClpSolverInterface * si =
1737                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
1738                ClpSimplex * clpSolver = si->getModelPtr();
1739                if (clpSolver->tightenPrimalBounds()!=0) {
1740                  std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
1741                  exit(1);
1742                }
1743                if (clpSolver->dualBound()==1.0e10) {
1744                  // user did not set - so modify
1745                  // get largest scaled away from bound
1746                  double largest=1.0e-12;
1747                  int numberRows = clpSolver->numberRows();
1748                  const double * rowPrimal = clpSolver->primalRowSolution();
1749                  const double * rowLower = clpSolver->rowLower();
1750                  const double * rowUpper = clpSolver->rowUpper();
1751                  const double * rowScale = clpSolver->rowScale();
1752                  int iRow;
1753                  for (iRow=0;iRow<numberRows;iRow++) {
1754                    double value = rowPrimal[iRow];
1755                    double above = value-rowLower[iRow];
1756                    double below = rowUpper[iRow]-value;
1757                    if (rowScale) {
1758                      double multiplier = rowScale[iRow];
1759                      above *= multiplier;
1760                      below *= multiplier;
1761                    }
1762                    if (above<1.0e12)
1763                      largest = CoinMax(largest,above);
1764                    if (below<1.0e12)
1765                      largest = CoinMax(largest,below);
1766                  }
1767                 
1768                  int numberColumns = clpSolver->numberColumns();
1769                  const double * columnPrimal = clpSolver->primalColumnSolution();
1770                  const double * columnLower = clpSolver->columnLower();
1771                  const double * columnUpper = clpSolver->columnUpper();
1772                  const double * columnScale = clpSolver->columnScale();
1773                  int iColumn;
1774                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1775                    double value = columnPrimal[iColumn];
1776                    double above = value-columnLower[iColumn];
1777                    double below = columnUpper[iColumn]-value;
1778                    if (columnScale) {
1779                      double multiplier = 1.0/columnScale[iColumn];
1780                      above *= multiplier;
1781                      below *= multiplier;
1782                    }
1783                    if (above<1.0e12)
1784                      largest = CoinMax(largest,above);
1785                    if (below<1.0e12)
1786                      largest = CoinMax(largest,below);
1787                  }
1788                  //if (!noPrinting)
1789                  //std::cout<<"Largest (scaled) away from bound "<<largest<<std::endl;
1790                  clpSolver->setDualBound(CoinMax(1.0001e8,CoinMin(1000.0*largest,1.00001e10)));
1791                }
1792                clpSolver->dual();  // clean up
1793              }
1794              // If user made settings then use them
1795              if (!defaultSettings) {
1796                OsiSolverInterface * solver = model.solver();
1797                if (!doScaling)
1798                  solver->setHintParam(OsiDoScale,false,OsiHintTry);
1799                OsiClpSolverInterface * si =
1800                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
1801                assert (si != NULL);
1802                // get clp itself
1803                ClpSimplex * modelC = si->getModelPtr();
1804                //if (modelC->tightenPrimalBounds()!=0) {
1805                //std::cout<<"Problem is infeasible!"<<std::endl;
1806                //break;
1807                //}
1808                // bounds based on continuous
1809                if (tightenFactor) {
1810                  if (modelC->tightenPrimalBounds(tightenFactor)!=0) {
1811                    std::cout<<"Problem is infeasible!"<<std::endl;
1812                    break;
1813                  }
1814                }
1815                if (djFix<1.0e20) {
1816                  // do some fixing
1817                  int numberColumns = modelC->numberColumns();
1818                  int i;
1819                  const char * type = modelC->integerInformation();
1820                  double * lower = modelC->columnLower();
1821                  double * upper = modelC->columnUpper();
1822                  double * solution = modelC->primalColumnSolution();
1823                  double * dj = modelC->dualColumnSolution();
1824                  int numberFixed=0;
1825                  for (i=0;i<numberColumns;i++) {
1826                    if (type[i]) {
1827                      double value = solution[i];
1828                      if (value<lower[i]+1.0e-5&&dj[i]>djFix) {
1829                        solution[i]=lower[i];
1830                        upper[i]=lower[i];
1831                        numberFixed++;
1832                      } else if (value>upper[i]-1.0e-5&&dj[i]<-djFix) {
1833                        solution[i]=upper[i];
1834                        lower[i]=upper[i];
1835                        numberFixed++;
1836                      }
1837                    }
1838                  }
1839                  printf("%d columns fixed\n",numberFixed);
1840                }
1841              }
1842              // See if we want preprocessing
1843              OsiSolverInterface * saveSolver=NULL;
1844              CglPreProcess process;
1845              delete babModel;
1846              babModel = new CbcModel(model);
1847              OsiSolverInterface * solver3 = clpSolver->clone();
1848              babModel->assignSolver(solver3);
1849              OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
1850              int numberChanged=0;
1851              if (clpSolver2->messageHandler()->logLevel())
1852                clpSolver2->messageHandler()->setLogLevel(1);
1853              if (logLevel>-1)
1854                clpSolver2->messageHandler()->setLogLevel(logLevel);
1855              lpSolver = clpSolver2->getModelPtr();
1856              if (lpSolver->factorizationFrequency()==200&&!miplib) {
1857                // User did not touch preset
1858                int numberRows = lpSolver->numberRows();
1859                const int cutoff1=10000;
1860                const int cutoff2=100000;
1861                const int base=75;
1862                const int freq0 = 50;
1863                const int freq1=200;
1864                const int freq2=400;
1865                const int maximum=1000;
1866                int frequency;
1867                if (numberRows<cutoff1)
1868                  frequency=base+numberRows/freq0;
1869                else if (numberRows<cutoff2)
1870                  frequency=base+cutoff1/freq0 + (numberRows-cutoff1)/freq1;
1871                else
1872                  frequency=base+cutoff1/freq0 + (cutoff2-cutoff1)/freq1 + (numberRows-cutoff2)/freq2;
1873                lpSolver->setFactorizationFrequency(CoinMin(maximum,frequency));
1874              }
1875              time2 = CoinCpuTime();
1876              totalTime += time2-time1;
1877              time1 = time2;
1878              double timeLeft = babModel->getMaximumSeconds();
1879              int numberOriginalColumns = babModel->solver()->getNumCols();
1880              if (preProcess==7) {
1881                // use strategy instead
1882                preProcess=0;
1883                useStrategy=true;
1884              }
1885              if (preProcess&&type==BAB) {
1886                // See if sos from mps file
1887                if (numberSOS==0&&clpSolver->numberSOS()&&doSOS) {
1888                  // SOS
1889                  numberSOS = clpSolver->numberSOS();
1890                  const CoinSet * setInfo = clpSolver->setInfo();
1891                  sosStart = new int [numberSOS+1];
1892                  sosType = new char [numberSOS];
1893                  int i;
1894                  int nTotal=0;
1895                  sosStart[0]=0;
1896                  for ( i=0;i<numberSOS;i++) {
1897                    int type = setInfo[i].setType();
1898                    int n=setInfo[i].numberEntries();
1899                    sosType[i]=type;
1900                    nTotal += n;
1901                    sosStart[i+1] = nTotal;
1902                  }
1903                  sosIndices = new int[nTotal];
1904                  sosReference = new double [nTotal];
1905                  for (i=0;i<numberSOS;i++) {
1906                    int n=setInfo[i].numberEntries();
1907                    const int * which = setInfo[i].which();
1908                    const double * weights = setInfo[i].weights();
1909                    int base = sosStart[i];
1910                    for (int j=0;j<n;j++) {
1911                      int k=which[j];
1912                      sosIndices[j+base]=k;
1913                      sosReference[j+base] = weights ? weights[j] : (double) j;
1914                    }
1915                  }
1916                }
1917                saveSolver=babModel->solver()->clone();
1918                /* Do not try and produce equality cliques and
1919                   do up to 10 passes */
1920                OsiSolverInterface * solver2;
1921                {
1922                  // Tell solver we are in Branch and Cut
1923                  saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;
1924                  // Default set of cut generators
1925                  CglProbing generator1;
1926                  generator1.setUsingObjective(true);
1927                  generator1.setMaxPass(3);
1928                  generator1.setMaxProbeRoot(saveSolver->getNumCols());
1929                  generator1.setMaxElements(100);
1930                  generator1.setMaxLookRoot(50);
1931                  generator1.setRowCuts(3);
1932                  // Add in generators
1933                  process.addCutGenerator(&generator1);
1934                  int translate[]={9999,0,0,-1,2,3,-2};
1935                  process.messageHandler()->setLogLevel(babModel->logLevel());
1936#ifdef COIN_HAS_ASL
1937                  if (info.numberSos&&doSOS&&usingAmpl) {
1938                    // SOS
1939                    numberSOS = info.numberSos;
1940                    sosStart = info.sosStart;
1941                    sosIndices = info.sosIndices;
1942                  }
1943#endif
1944                  if (numberSOS&&doSOS) {
1945                    // SOS
1946                    int numberColumns = saveSolver->getNumCols();
1947                    char * prohibited = new char[numberColumns];
1948                    memset(prohibited,0,numberColumns);
1949                    int n=sosStart[numberSOS];
1950                    for (int i=0;i<n;i++) {
1951                      int iColumn = sosIndices[i];
1952                      prohibited[iColumn]=1;
1953                    }
1954                    process.passInProhibited(prohibited,numberColumns);
1955                    delete [] prohibited;
1956                  }
1957                  solver2 = process.preProcessNonDefault(*saveSolver,translate[preProcess],10,
1958                                                         tunePreProcess);
1959                  // Tell solver we are not in Branch and Cut
1960                  saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
1961                  if (solver2)
1962                    solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
1963                }
1964#ifdef COIN_HAS_ASL
1965                if (!solver2&&usingAmpl) {
1966                  // infeasible
1967                  info.problemStatus=1;
1968                  info.objValue = 1.0e100;
1969                  sprintf(info.buffer,"infeasible/unbounded by pre-processing");
1970                  info.primalSolution=NULL;
1971                  info.dualSolution=NULL;
1972                  break;
1973                }
1974#endif
1975                if (!noPrinting) {
1976                  if (!solver2) {
1977                    printf("Pre-processing says infeasible or unbounded\n");
1978                    break;
1979                  } else {
1980                    printf("processed model has %d rows, %d columns and %d elements\n",
1981                           solver2->getNumRows(),solver2->getNumCols(),solver2->getNumElements());
1982                  }
1983                }
1984                //solver2->resolve();
1985                if (preProcess==2) {
1986                  OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (solver2);
1987                  ClpSimplex * lpSolver = clpSolver2->getModelPtr();
1988                  lpSolver->writeMps("presolved.mps",0,1,lpSolver->optimizationDirection());
1989                  printf("Preprocessed model (minimization) on presolved.mps\n");
1990                }
1991                // we have to keep solver2 so pass clone
1992                solver2 = solver2->clone();
1993                babModel->assignSolver(solver2);
1994                babModel->initialSolve();
1995                babModel->setMaximumSeconds(timeLeft-(CoinCpuTime()-time1));
1996              }
1997              // now tighten bounds
1998              if (!miplib) {
1999                OsiClpSolverInterface * si =
2000                  dynamic_cast<OsiClpSolverInterface *>(babModel->solver()) ;
2001                assert (si != NULL);
2002                // get clp itself
2003                ClpSimplex * modelC = si->getModelPtr();
2004                if (noPrinting)
2005                  modelC->setLogLevel(0);
2006                if (modelC->tightenPrimalBounds()!=0) {
2007                  std::cout<<"Problem is infeasible!"<<std::endl;
2008                  break;
2009                }
2010                modelC->dual();
2011              }
2012              if (debugValues) {
2013                // for debug
2014                std::string problemName ;
2015                babModel->solver()->getStrParam(OsiProbName,problemName) ;
2016                babModel->solver()->activateRowCutDebugger(problemName.c_str()) ;
2017                twomirGen.probname_=strdup(problemName.c_str());
2018                // checking seems odd
2019                //redsplitGen.set_given_optsol(babModel->solver()->getRowCutDebuggerAlways()->optimalSolution(),
2020                //                         babModel->getNumCols());
2021              }
2022              if (useCosts) {
2023                int numberColumns = babModel->getNumCols();
2024                int * sort = new int[numberColumns];
2025                double * dsort = new double[numberColumns];
2026                int * priority = new int [numberColumns];
2027                const double * objective = babModel->getObjCoefficients();
2028                const double * lower = babModel->getColLower() ;
2029                const double * upper = babModel->getColUpper() ;
2030                const CoinPackedMatrix * matrix = babModel->solver()->getMatrixByCol();
2031                const int * columnLength = matrix->getVectorLengths();
2032                int iColumn;
2033                int n=0;
2034                for (iColumn=0;iColumn<numberColumns;iColumn++) {
2035                  if (babModel->isInteger(iColumn)) {
2036                    sort[n]=n;
2037                    if (useCosts==1)
2038                      dsort[n++]=-fabs(objective[iColumn]);
2039                    else if (useCosts==2)
2040                      dsort[n++]=iColumn;
2041                    else if (useCosts==3)
2042                      dsort[n++]=upper[iColumn]-lower[iColumn];
2043                    else if (useCosts==4)
2044                      dsort[n++]=-(upper[iColumn]-lower[iColumn]);
2045                    else if (useCosts==5)
2046                      dsort[n++]=-columnLength[iColumn];
2047                  }
2048                }
2049                CoinSort_2(dsort,dsort+n,sort);
2050                int level=0;
2051                double last = -1.0e100;
2052                for (int i=0;i<n;i++) {
2053                  int iPut=sort[i];
2054                  if (dsort[i]!=last) {
2055                    level++;
2056                    last=dsort[i];
2057                  }
2058                  priority[iPut]=level;
2059                }
2060                babModel->passInPriorities( priority,false);
2061                delete [] priority;
2062                delete [] sort;
2063                delete [] dsort;
2064              }
2065              // FPump done first as it only works if no solution
2066              CbcHeuristicFPump heuristic4(*babModel);
2067              if (useFpump) {
2068                heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
2069                int pumpTune=parameters[whichParam(FPUMPTUNE,numberParameters,parameters)].intValue();
2070                if (pumpTune>0) {
2071                  /*
2072                    >=10000000 for using obj
2073                    >=1000000 use as accumulate switch
2074                    >=1000 use index+2 as number of large loops
2075                    >=100 use 0.05 objvalue as increment
2076                    >=10 use +0.1 objvalue for cutoff (add)
2077                    1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
2078                    4 and static continuous, 5 as 3 but no internal integers
2079                  */
2080                  double value = babModel->solver()->getObjSense()*babModel->solver()->getObjValue();
2081                  int w = pumpTune/10;
2082                  int c = w % 10;
2083                  w /= 10;
2084                  int i = w % 10;
2085                  w /= 10;
2086                  int r = w;
2087                  int accumulate = r/1000;
2088                  r -= 1000*accumulate;
2089                  if (accumulate>=10) {
2090                    int which = accumulate/10;
2091                    accumulate -= 10*which;
2092                    which--;
2093                    // weights and factors
2094                    double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
2095                    double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
2096                    heuristic4.setInitialWeight(weight[which]);
2097                    heuristic4.setWeightFactor(factor[which]);
2098                  }
2099                  // fake cutoff
2100                  printf("Setting ");
2101                  if (c) {
2102                    double cutoff;
2103                    babModel->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
2104                    cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
2105                    heuristic4.setFakeCutoff(cutoff);
2106                    printf("fake cutoff of %g ",cutoff);
2107                  }
2108                  if (i||r) {
2109                    // also set increment
2110                    heuristic4.setAbsoluteIncrement((0.01*i+0.005)*(fabs(value)+1.0e-12));
2111                    heuristic4.setAccumulate(accumulate);
2112                    heuristic4.setMaximumRetries(r+2);
2113                    if (i)
2114                      printf("increment of %g ",heuristic4.absoluteIncrement());
2115                    if (accumulate)
2116                      printf("accumulate of %d ",accumulate);
2117                    printf("%d retries ",r+2);
2118                  }
2119                  pumpTune = pumpTune%100;
2120                  printf("and setting when to %d\n",pumpTune+10);
2121                  heuristic4.setWhen(pumpTune+10);
2122                }
2123                babModel->addHeuristic(&heuristic4);
2124              }
2125              if (!miplib) {
2126                CbcRounding heuristic1(*babModel);
2127                if (useRounding)
2128                  babModel->addHeuristic(&heuristic1) ;
2129                CbcHeuristicLocal heuristic2(*babModel);
2130                heuristic2.setSearchType(1);
2131                if (useCombine)
2132                  babModel->addHeuristic(&heuristic2);
2133                CbcHeuristicGreedyCover heuristic3(*babModel);
2134                CbcHeuristicGreedyEquality heuristic3a(*babModel);
2135                if (useGreedy) {
2136                  babModel->addHeuristic(&heuristic3);
2137                  babModel->addHeuristic(&heuristic3a);
2138                }
2139                if (useLocalTree) {
2140                  CbcTreeLocal localTree(babModel,NULL,10,0,0,10000,2000);
2141                  babModel->passInTreeHandler(localTree);
2142                }
2143              }
2144              if (type==MIPLIB) {
2145                if (babModel->numberStrong()==5&&babModel->numberBeforeTrust()==5) 
2146                  babModel->setNumberBeforeTrust(50);
2147              }
2148              // add cut generators if wanted
2149              int switches[20];
2150              int numberGenerators=0;
2151              int translate[]={-100,-1,-99,-98,1,1,1,1};
2152              if (probingAction) {
2153                if (probingAction==5||probingAction==7)
2154                  probingGen.setRowCuts(-3); // strengthening etc just at root
2155                if (probingAction==6||probingAction==7) {
2156                  // Number of unsatisfied variables to look at
2157                  probingGen.setMaxProbe(1000);
2158                  probingGen.setMaxProbeRoot(1000);
2159                  // How far to follow the consequences
2160                  probingGen.setMaxLook(50);
2161                  probingGen.setMaxLookRoot(50);
2162                }
2163                babModel->addCutGenerator(&probingGen,translate[probingAction],"Probing");
2164                switches[numberGenerators++]=0;
2165              }
2166              if (gomoryAction) {
2167                babModel->addCutGenerator(&gomoryGen,translate[gomoryAction],"Gomory");
2168                switches[numberGenerators++]=-1;
2169              }
2170              if (knapsackAction) {
2171                babModel->addCutGenerator(&knapsackGen,translate[knapsackAction],"Knapsack");
2172                switches[numberGenerators++]=0;
2173              }
2174              if (redsplitAction) {
2175                babModel->addCutGenerator(&redsplitGen,translate[redsplitAction],"Reduce-and-split");
2176                switches[numberGenerators++]=1;
2177              }
2178              if (cliqueAction) {
2179                babModel->addCutGenerator(&cliqueGen,translate[cliqueAction],"Clique");
2180                switches[numberGenerators++]=0;
2181              }
2182              if (mixedAction) {
2183                babModel->addCutGenerator(&mixedGen,translate[mixedAction],"MixedIntegerRounding2");
2184                switches[numberGenerators++]=-1;
2185              }
2186              if (flowAction) {
2187                babModel->addCutGenerator(&flowGen,translate[flowAction],"FlowCover");
2188                switches[numberGenerators++]=1;
2189              }
2190              if (twomirAction) {
2191                babModel->addCutGenerator(&twomirGen,translate[twomirAction],"TwoMirCuts");
2192                switches[numberGenerators++]=1;
2193              }
2194              if (landpAction) {
2195                babModel->addCutGenerator(&landpGen,translate[landpAction],"LiftAndProject");
2196                switches[numberGenerators++]=1;
2197              }
2198              if (storedCuts) 
2199                babModel->setSpecialOptions(babModel->specialOptions()|64);
2200              // Say we want timings
2201              numberGenerators = babModel->numberCutGenerators();
2202              int iGenerator;
2203              int cutDepth=
2204                parameters[whichParam(CUTDEPTH,numberParameters,parameters)].intValue();
2205              for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
2206                CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
2207                int howOften = generator->howOften();
2208                if (howOften==-98||howOften==-99) 
2209                  generator->setSwitchOffIfLessThan(switches[iGenerator]);
2210                generator->setTiming(true);
2211                if (cutDepth>=0)
2212                  generator->setWhatDepth(cutDepth) ;
2213              }
2214              // Could tune more
2215              if (!miplib) {
2216                babModel->setMinimumDrop(min(5.0e-2,
2217                                             fabs(babModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
2218                if (cutPass==-1234567) {
2219                  if (babModel->getNumCols()<500)
2220                    babModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
2221                  else if (babModel->getNumCols()<5000)
2222                    babModel->setMaximumCutPassesAtRoot(100); // use minimum drop
2223                  else
2224                    babModel->setMaximumCutPassesAtRoot(20);
2225                } else {
2226                  babModel->setMaximumCutPassesAtRoot(cutPass);
2227                }
2228                babModel->setMaximumCutPasses(1);
2229              }
2230              // Do more strong branching if small
2231              //if (babModel->getNumCols()<5000)
2232              //babModel->setNumberStrong(20);
2233              // Switch off strong branching if wanted
2234              //if (babModel->getNumCols()>10*babModel->getNumRows())
2235              //babModel->setNumberStrong(0);
2236              if (!noPrinting) {
2237                int iLevel = parameters[log].intValue();
2238                if (iLevel<0) {
2239                  babModel->setPrintingMode(1);
2240                  iLevel = -iLevel;
2241                }
2242                babModel->messageHandler()->setLogLevel(iLevel);
2243                if (babModel->getNumCols()>2000||babModel->getNumRows()>1500||
2244                    babModel->messageHandler()->logLevel()>1)
2245                  babModel->setPrintFrequency(100);
2246              }
2247             
2248              babModel->solver()->setIntParam(OsiMaxNumIterationHotStart,
2249                    parameters[whichParam(MAXHOTITS,numberParameters,parameters)].intValue());
2250              OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
2251              // go faster stripes
2252              if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
2253                osiclp->setupForRepeatedUse(2,parameters[slog].intValue());
2254              } else {
2255                osiclp->setupForRepeatedUse(0,parameters[slog].intValue());
2256              }
2257              double increment=babModel->getCutoffIncrement();;
2258              int * changed = NULL;
2259              if (!miplib)
2260                changed=analyze( osiclp,numberChanged,increment,false);
2261              if (debugValues) {
2262                if (numberDebugValues==babModel->getNumCols()) {
2263                  // for debug
2264                  babModel->solver()->activateRowCutDebugger(debugValues) ;
2265                } else {
2266                  printf("debug file has incorrect number of columns\n");
2267                }
2268              }
2269              babModel->setCutoffIncrement(CoinMax(babModel->getCutoffIncrement(),increment));
2270              // Turn this off if you get problems
2271              // Used to be automatically set
2272              int mipOptions = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue();
2273              if (mipOptions!=(128|64|1))
2274                printf("mip options %d\n",mipOptions);
2275              osiclp->setSpecialOptions(mipOptions);
2276              if (gapRatio < 1.0e100) {
2277                double value = babModel->solver()->getObjValue() ;
2278                double value2 = gapRatio*(1.0e-5+fabs(value)) ;
2279                babModel->setAllowableGap(value2) ;
2280                std::cout << "Continuous " << value
2281                          << ", so allowable gap set to "
2282                          << value2 << std::endl ;
2283              }
2284              // probably faster to use a basis to get integer solutions
2285              babModel->setSpecialOptions(babModel->specialOptions()|2);
2286              currentBranchModel = babModel;
2287              OsiSolverInterface * strengthenedModel=NULL;
2288              if (type==BAB||type==MIPLIB) {
2289                int moreMipOptions = parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].intValue();
2290                if (moreMipOptions>=0) {
2291                  printf("more mip options %d\n",moreMipOptions);
2292                  if (((moreMipOptions+1)%1000000)!=0)
2293                    babModel->setSearchStrategy(moreMipOptions%1000000);
2294                  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
2295                  // go faster stripes
2296                  if( moreMipOptions >=999999) {
2297                    if (osiclp) {
2298                      int save = osiclp->specialOptions();
2299                      osiclp->setupForRepeatedUse(2,0);
2300                      osiclp->setSpecialOptions(save|osiclp->specialOptions());
2301                    }
2302                  } 
2303                }
2304              }
2305              if (type==BAB) {
2306#ifdef COIN_HAS_ASL
2307                if (usingAmpl) {
2308                  priorities=info.priorities;
2309                  branchDirection=info.branchDirection;
2310                  pseudoDown=info.pseudoDown;
2311                  pseudoUp=info.pseudoUp;
2312                  solutionIn=info.primalSolution;
2313                  prioritiesIn = info.priorities;
2314                  if (info.numberSos&&doSOS) {
2315                    // SOS
2316                    numberSOS = info.numberSos;
2317                    sosStart = info.sosStart;
2318                    sosIndices = info.sosIndices;
2319                    sosType = info.sosType;
2320                    sosReference = info.sosReference;
2321                    sosPriority = info.sosPriority;
2322                  }
2323                }
2324#endif               
2325                const int * originalColumns = preProcess ? process.originalColumns() : NULL;
2326                if (solutionIn&&useSolution) {
2327                  if (preProcess) {
2328                    int numberColumns = babModel->getNumCols();
2329                    // extend arrays in case SOS
2330                    int n = originalColumns[numberColumns-1]+1;
2331                    int nSmaller = CoinMin(n,numberOriginalColumns);
2332                    double * solutionIn2 = new double [n];
2333                    int * prioritiesIn2 = new int[n];
2334                    int i;
2335                    for (i=0;i<nSmaller;i++) {
2336                      solutionIn2[i]=solutionIn[i];
2337                      prioritiesIn2[i]=prioritiesIn[i];
2338                    }
2339                    for (;i<n;i++) {
2340                      solutionIn2[i]=0.0;
2341                      prioritiesIn2[i]=1000000;
2342                    }
2343                    int iLast=-1;
2344                    for (i=0;i<numberColumns;i++) {
2345                      int iColumn = originalColumns[i];
2346                      assert (iColumn>iLast);
2347                      iLast=iColumn;
2348                      solutionIn2[i]=solutionIn2[iColumn];
2349                      if (prioritiesIn)
2350                        prioritiesIn2[i]=prioritiesIn2[iColumn];
2351                    }
2352                    babModel->setHotstartSolution(solutionIn2,prioritiesIn2);
2353                    delete [] solutionIn2;
2354                    delete [] prioritiesIn2;
2355                  } else {
2356                    babModel->setHotstartSolution(solutionIn,prioritiesIn);
2357                  }
2358                }
2359                if (preProcess&&process.numberSOS()) {
2360                  int numberSOS = process.numberSOS();
2361                  int numberIntegers = babModel->numberIntegers();
2362                  /* model may not have created objects
2363                     If none then create
2364                  */
2365                  if (!numberIntegers||!babModel->numberObjects()) {
2366                    int type = (pseudoUp) ? 1 : 0;
2367                    babModel->findIntegers(true,type);
2368                    numberIntegers = babModel->numberIntegers();
2369                  }
2370                  OsiObject ** oldObjects = babModel->objects();
2371                  // Do sets and priorities
2372                  OsiObject ** objects = new OsiObject * [numberSOS];
2373                  // set old objects to have low priority
2374                  int numberOldObjects = babModel->numberObjects();
2375                  int numberColumns = babModel->getNumCols();
2376                  for (int iObj = 0;iObj<numberOldObjects;iObj++) {
2377                    oldObjects[iObj]->setPriority(numberColumns+1);
2378                    int iColumn = oldObjects[iObj]->columnNumber();
2379                    assert (iColumn>=0);
2380                    if (iColumn>=numberOriginalColumns)
2381                      continue;
2382                    if (originalColumns)
2383                      iColumn = originalColumns[iColumn];
2384                    if (branchDirection) {
2385                      CbcSimpleInteger * obj =
2386                        dynamic_cast <CbcSimpleInteger *>(oldObjects[iObj]) ;
2387                      if (obj) { 
2388                        obj->setPreferredWay(branchDirection[iColumn]);
2389                      } else {
2390                        CbcObject * obj =
2391                          dynamic_cast <CbcObject *>(oldObjects[iObj]) ;
2392                        assert (obj);
2393                        obj->setPreferredWay(branchDirection[iColumn]);
2394                      }
2395                    }
2396                    if (pseudoUp) {
2397                      CbcSimpleIntegerPseudoCost * obj1a =
2398                        dynamic_cast <CbcSimpleIntegerPseudoCost *>(oldObjects[iObj]) ;
2399                      assert (obj1a);
2400                      if (pseudoDown[iColumn]>0.0)
2401                        obj1a->setDownPseudoCost(pseudoDown[iColumn]);
2402                      if (pseudoUp[iColumn]>0.0)
2403                        obj1a->setUpPseudoCost(pseudoUp[iColumn]);
2404                    }
2405                  }
2406                  const int * starts = process.startSOS();
2407                  const int * which = process.whichSOS();
2408                  const int * type = process.typeSOS();
2409                  const double * weight = process.weightSOS();
2410                  int iSOS;
2411                  for (iSOS =0;iSOS<numberSOS;iSOS++) {
2412                    int iStart = starts[iSOS];
2413                    int n=starts[iSOS+1]-iStart;
2414                    objects[iSOS] = new CbcSOS(babModel,n,which+iStart,weight+iStart,
2415                                               iSOS,type[iSOS]);
2416                    // branch on long sets first
2417                    objects[iSOS]->setPriority(numberColumns-n);
2418                  }
2419                  babModel->addObjects(numberSOS,objects);
2420                  for (iSOS=0;iSOS<numberSOS;iSOS++)
2421                    delete objects[iSOS];
2422                  delete [] objects;
2423                } else if (priorities||branchDirection||pseudoDown||pseudoUp||numberSOS) {
2424                  // do anyway for priorities etc
2425                  int numberIntegers = babModel->numberIntegers();
2426                  /* model may not have created objects
2427                     If none then create
2428                  */
2429                  if (!numberIntegers||!babModel->numberObjects()) {
2430                    int type = (pseudoUp) ? 1 : 0;
2431                    babModel->findIntegers(true,type);
2432                  }
2433                  if (numberSOS) {
2434                    // Do sets and priorities
2435                    OsiObject ** objects = new OsiObject * [numberSOS];
2436                    int iSOS;
2437                    if (originalColumns) {
2438                      // redo sequence numbers
2439                      int numberColumns = babModel->getNumCols();
2440                      int nOld = originalColumns[numberColumns-1]+1;
2441                      int * back = new int[nOld];
2442                      int i;
2443                      for (i=0;i<nOld;i++)
2444                        back[i]=-1;
2445                      for (i=0;i<numberColumns;i++)
2446                        back[originalColumns[i]]=i;
2447                      // Really need better checks
2448                      int nMissing=0;
2449                      int n=sosStart[numberSOS];
2450                      for (i=0;i<n;i++) {
2451                        int iColumn = sosIndices[i];
2452                        int jColumn = back[iColumn];
2453                        if (jColumn>=0) 
2454                          sosIndices[i] = jColumn;
2455                        else 
2456                          nMissing++;
2457                      }
2458                      delete [] back;
2459                      if (nMissing)
2460                        printf("%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
2461                    }
2462                    for (iSOS =0;iSOS<numberSOS;iSOS++) {
2463                      int iStart = sosStart[iSOS];
2464                      int n=sosStart[iSOS+1]-iStart;
2465                      objects[iSOS] = new CbcSOS(babModel,n,sosIndices+iStart,sosReference+iStart,
2466                                                 iSOS,sosType[iSOS]);
2467                      if (sosPriority)
2468                        objects[iSOS]->setPriority(sosPriority[iSOS]);
2469                      else if (!prioritiesIn)
2470                        objects[iSOS]->setPriority(10);  // rather than 1000
2471                    }
2472                    // delete any existing SOS objects
2473                    int numberObjects=babModel->numberObjects();
2474                    OsiObject ** oldObjects=babModel->objects();
2475                    int nNew=0;
2476                    for (int i=0;i<numberObjects;i++) {
2477                      OsiObject * objThis = oldObjects[i];
2478                      CbcSOS * obj1 =
2479                        dynamic_cast <CbcSOS *>(objThis) ;
2480                      OsiSOS * obj2 =
2481                        dynamic_cast <OsiSOS *>(objThis) ;
2482                      if (!obj1&&!obj2) {
2483                        oldObjects[nNew++]=objThis;
2484                      } else {
2485                        delete objThis;
2486                      }
2487                    }
2488                    babModel->setNumberObjects(nNew);
2489                    babModel->addObjects(numberSOS,objects);
2490                    for (iSOS=0;iSOS<numberSOS;iSOS++)
2491                      delete objects[iSOS];
2492                    delete [] objects;
2493                  }
2494                  OsiObject ** objects = babModel->objects();
2495                  int numberObjects = babModel->numberObjects();
2496                  for (int iObj = 0;iObj<numberObjects;iObj++) {
2497                    // skip sos
2498                    CbcSOS * objSOS =
2499                      dynamic_cast <CbcSOS *>(objects[iObj]) ;
2500                    if (objSOS)
2501                      continue;
2502                    int iColumn = objects[iObj]->columnNumber();
2503                    assert (iColumn>=0);
2504                    if (originalColumns)
2505                      iColumn = originalColumns[iColumn];
2506                    if (branchDirection) {
2507                      CbcSimpleInteger * obj =
2508                        dynamic_cast <CbcSimpleInteger *>(objects[iObj]) ;
2509                      if (obj) { 
2510                        obj->setPreferredWay(branchDirection[iColumn]);
2511                      } else {
2512                        CbcObject * obj =
2513                          dynamic_cast <CbcObject *>(objects[iObj]) ;
2514                        assert (obj);
2515                        obj->setPreferredWay(branchDirection[iColumn]);
2516                      }
2517                    }
2518                    if (priorities) {
2519                      int iPriority = priorities[iColumn];
2520                      if (iPriority>0)
2521                        objects[iObj]->setPriority(iPriority);
2522                    }
2523                    if (pseudoUp&&pseudoUp[iColumn]) {
2524                      CbcSimpleIntegerPseudoCost * obj1a =
2525                        dynamic_cast <CbcSimpleIntegerPseudoCost *>(objects[iObj]) ;
2526                      assert (obj1a);
2527                      if (pseudoDown[iColumn]>0.0)
2528                        obj1a->setDownPseudoCost(pseudoDown[iColumn]);
2529                      if (pseudoUp[iColumn]>0.0)
2530                        obj1a->setUpPseudoCost(pseudoUp[iColumn]);
2531                    }
2532                  }
2533                }
2534                int statistics = (printOptions>0) ? printOptions: 0;
2535#ifdef COIN_HAS_ASL
2536                if (!usingAmpl) {
2537#endif
2538                  free(priorities);
2539                  priorities=NULL;
2540                  free(branchDirection);
2541                  branchDirection=NULL;
2542                  free(pseudoDown);
2543                  pseudoDown=NULL;
2544                  free(pseudoUp);
2545                  pseudoUp=NULL;
2546                  free(solutionIn);
2547                  solutionIn=NULL;
2548                  free(prioritiesIn);
2549                  prioritiesIn=NULL;
2550                  free(sosStart);
2551                  sosStart=NULL;
2552                  free(sosIndices);
2553                  sosIndices=NULL;
2554                  free(sosType);
2555                  sosType=NULL;
2556                  free(sosReference);
2557                  sosReference=NULL;
2558                  free(sosPriority);
2559                  sosPriority=NULL;
2560#ifdef COIN_HAS_ASL
2561                }
2562#endif               
2563                if (nodeStrategy) {
2564                  // change default
2565                  if (nodeStrategy>1) {
2566                    // up or down
2567                    int way = ((nodeStrategy%1)==1) ? -1 : +1;
2568                    babModel->setPreferredWay(way);
2569#if 0
2570                    OsiObject ** objects = babModel->objects();
2571                    int numberObjects = babModel->numberObjects();
2572                    for (int iObj = 0;iObj<numberObjects;iObj++) {
2573                      CbcObject * obj =
2574                        dynamic_cast <CbcObject *>(objects[iObj]) ;
2575                      assert (obj);
2576                      obj->setPreferredWay(way);
2577                    }
2578#endif
2579                  }
2580                  if (nodeStrategy==1||nodeStrategy>3) {
2581                    // depth
2582                    CbcCompareDefault compare;
2583                    compare.setWeight(-3.0);
2584                    babModel->setNodeComparison(compare);
2585                  }
2586                }
2587                if (cppValue>=0) {
2588                  int prepro = useStrategy ? -1 : preProcess;
2589                  // generate code
2590                  FILE * fp = fopen("user_driver.cpp","w");
2591                  if (fp) {
2592                    // generate enough to do BAB
2593                    babModel->generateCpp(fp,1);
2594                    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
2595                    // Make general so do factorization
2596                    int factor = osiclp->getModelPtr()->factorizationFrequency();
2597                    osiclp->getModelPtr()->setFactorizationFrequency(200);
2598                    osiclp->generateCpp(fp);
2599                    osiclp->getModelPtr()->setFactorizationFrequency(factor);
2600                    //solveOptions.generateCpp(fp);
2601                    fclose(fp);
2602                    // now call generate code
2603                    generateCode(babModel,"user_driver.cpp",cppValue,prepro);
2604                  } else {
2605                    std::cout<<"Unable to open file user_driver.cpp"<<std::endl;
2606                  }
2607                }
2608                if (useStrategy) {
2609                  CbcStrategyDefault strategy(true,babModel->numberStrong(),babModel->numberBeforeTrust());
2610                  strategy.setupPreProcessing(1);
2611                  babModel->setStrategy(strategy);
2612                }
2613                int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
2614                if (testOsiOptions>=0) {
2615                  printf("Testing OsiObject options %d\n",testOsiOptions);
2616                  CbcBranchDefaultDecision decision;
2617                  babModel->solver()->findIntegersAndSOS(false);
2618                  //OsiChooseVariable choose(babModel->solver());
2619                  OsiChooseStrong choose(babModel->solver());
2620                  choose.setNumberBeforeTrusted(babModel->numberBeforeTrust());
2621                  choose.setNumberStrong(babModel->numberStrong());
2622                  choose.setShadowPriceMode(testOsiOptions);
2623                  decision.setChooseMethod(choose);
2624                  babModel->setBranchingMethod(decision);
2625                }
2626                checkSOS(babModel, babModel->solver());
2627                babModel->branchAndBound(statistics);
2628                checkSOS(babModel, babModel->solver());
2629              } else if (type==MIPLIB) {
2630                CbcStrategyDefault strategy(true,babModel->numberStrong(),babModel->numberBeforeTrust());
2631                // Set up pre-processing
2632                int translate2[]={9999,1,1,3,2,4,5};
2633                if (preProcess)
2634                  strategy.setupPreProcessing(translate2[preProcess]);
2635                babModel->setStrategy(strategy);
2636                CbcClpUnitTest(*babModel);
2637                goodModel=false;
2638                break;
2639              } else {
2640                strengthenedModel = babModel->strengthenedModel();
2641              }
2642              currentBranchModel = NULL;
2643              osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
2644              if (debugFile=="createAfterPre"&&babModel->bestSolution()) {
2645                lpSolver = osiclp->getModelPtr();
2646                //move best solution (should be there -- but ..)
2647                int n = lpSolver->getNumCols();
2648                memcpy(lpSolver->primalColumnSolution(),babModel->bestSolution(),n*sizeof(double));
2649                saveSolution(osiclp->getModelPtr(),"debug.file");
2650              }
2651              if (!noPrinting) {
2652                // Print more statistics
2653                std::cout<<"Cuts at root node changed objective from "<<babModel->getContinuousObjective()
2654                         <<" to "<<babModel->rootObjectiveAfterCuts()<<std::endl;
2655               
2656                numberGenerators = babModel->numberCutGenerators();
2657                for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
2658                  CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
2659                  std::cout<<generator->cutGeneratorName()<<" was tried "
2660                           <<generator->numberTimesEntered()<<" times and created "
2661                           <<generator->numberCutsInTotal()<<" cuts of which "
2662                           <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
2663                  if (generator->timing())
2664                    std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
2665                  else
2666                    std::cout<<std::endl;
2667                }
2668              }
2669              time2 = CoinCpuTime();
2670              totalTime += time2-time1;
2671              // For best solution
2672              double * bestSolution = NULL;
2673              if (babModel->getMinimizationObjValue()<1.0e50&&type==BAB) {
2674                // post process
2675                if (preProcess) {
2676                  int n = saveSolver->getNumCols();
2677                  bestSolution = new double [n];
2678                  process.postProcess(*babModel->solver());
2679                  // Solution now back in saveSolver
2680                  babModel->assignSolver(saveSolver);
2681                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
2682                } else {
2683                  int n = babModel->solver()->getNumCols();
2684                  bestSolution = new double [n];
2685                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
2686                }
2687                checkSOS(babModel, babModel->solver());
2688              }
2689              if (type==STRENGTHEN&&strengthenedModel)
2690                clpSolver = dynamic_cast< OsiClpSolverInterface*> (strengthenedModel);
2691              lpSolver = clpSolver->getModelPtr();
2692              if (numberChanged) {
2693                for (int i=0;i<numberChanged;i++) {
2694                  int iColumn=changed[i];
2695                  clpSolver->setContinuous(iColumn);
2696                }
2697                delete [] changed;
2698              }
2699              if (type==BAB) {
2700                //move best solution (should be there -- but ..)
2701                int n = lpSolver->getNumCols();
2702                if (bestSolution)
2703                  memcpy(lpSolver->primalColumnSolution(),bestSolution,n*sizeof(double));
2704                if (debugFile=="create"&&bestSolution) {
2705                  saveSolution(lpSolver,"debug.file");
2706                }
2707                delete [] bestSolution;
2708                std::string statusName[]={"Finished","Stopped on ","Difficulties",
2709                                          "","","User ctrl-c"};
2710                std::string minor[]={"","","gap","nodes","time","","solutions","user ctrl-c"};
2711                int iStat = babModel->status();
2712                int iStat2 = babModel->secondaryStatus();
2713                if (!noPrinting)
2714                  std::cout<<"Result - "<<statusName[iStat]<<minor[iStat2]
2715                           <<" objective "<<babModel->getObjValue()<<
2716                    " after "<<babModel->getNodeCount()<<" nodes and "
2717                           <<babModel->getIterationCount()<<
2718                    " iterations - took "<<time2-time1<<" seconds"<<std::endl;
2719#ifdef COIN_HAS_ASL
2720                if (usingAmpl) {
2721                  double value = babModel->getObjValue()*lpSolver->getObjSense();
2722                  char buf[300];
2723                  int pos=0;
2724                  if (iStat==0) {
2725                    if (babModel->getObjValue()<1.0e40) {
2726                      pos += sprintf(buf+pos,"optimal," );
2727                    } else {
2728                      // infeasible
2729                      iStat=1;
2730                      pos += sprintf(buf+pos,"infeasible,");
2731                    }
2732                  } else if (iStat==1) {
2733                    if (iStat2!=6)
2734                      iStat=3;
2735                    else
2736                      iStat=4;
2737                    pos += sprintf(buf+pos,"stopped on %s,",minor[iStat2].c_str());
2738                  } else if (iStat==2) {
2739                    iStat = 7;
2740                    pos += sprintf(buf+pos,"stopped on difficulties,");
2741                  } else if (iStat==5) {
2742                    iStat = 3;
2743                    pos += sprintf(buf+pos,"stopped on ctrl-c,");
2744                  } else {
2745                    pos += sprintf(buf+pos,"status unknown,");
2746                    iStat=6;
2747                  }
2748                  info.problemStatus=iStat;
2749                  info.objValue = value;
2750                  if (babModel->getObjValue()<1.0e40) 
2751                    pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
2752                                   value);
2753                  sprintf(buf+pos,"\n%d nodes, %d iterations",
2754                          babModel->getNodeCount(),
2755                          babModel->getIterationCount());
2756                  if (bestSolution) {
2757                    free(info.primalSolution);
2758                    info.primalSolution = (double *) malloc(n*sizeof(double));
2759                    CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution);
2760                    int numberRows = lpSolver->numberRows();
2761                    free(info.dualSolution);
2762                    info.dualSolution = (double *) malloc(numberRows*sizeof(double));
2763                    CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution);
2764                  } else {
2765                    info.primalSolution=NULL;
2766                    info.dualSolution=NULL;
2767                  }
2768                  // put buffer into info
2769                  strcpy(info.buffer,buf);
2770                }
2771#endif
2772              } else {
2773                std::cout<<"Model strengthened - now has "<<clpSolver->getNumRows()
2774                         <<" rows"<<std::endl;
2775              }
2776              time1 = time2;
2777              delete babModel;
2778              babModel=NULL;
2779            } else {
2780              std::cout << "** Current model not valid" << std::endl ; 
2781            }
2782            break ;
2783          case IMPORT:
2784            {
2785#ifdef COIN_HAS_ASL
2786              if (!usingAmpl) {
2787#endif
2788                free(priorities);
2789                priorities=NULL;
2790                free(branchDirection);
2791                branchDirection=NULL;
2792                free(pseudoDown);
2793                pseudoDown=NULL;
2794                free(pseudoUp);
2795                pseudoUp=NULL;
2796                free(solutionIn);
2797                solutionIn=NULL;
2798                free(prioritiesIn);
2799                prioritiesIn=NULL;
2800                free(sosStart);
2801                sosStart=NULL;
2802                free(sosIndices);
2803                sosIndices=NULL;
2804                free(sosType);
2805                sosType=NULL;
2806                free(sosReference);
2807                sosReference=NULL;
2808                free(sosPriority);
2809                sosPriority=NULL;
2810#ifdef COIN_HAS_ASL
2811              }
2812#endif               
2813              delete babModel;
2814              babModel=NULL;
2815              // get next field
2816              field = CoinReadGetString(argc,argv);
2817              if (field=="$") {
2818                field = parameters[iParam].stringValue();
2819              } else if (field=="EOL") {
2820                parameters[iParam].printString();
2821                break;
2822              } else {
2823                parameters[iParam].setStringValue(field);
2824              }
2825              std::string fileName;
2826              bool canOpen=false;
2827              if (field=="-") {
2828                // stdin
2829                canOpen=true;
2830                fileName = "-";
2831              } else {
2832                bool absolutePath;
2833                if (dirsep=='/') {
2834                  // non Windows (or cygwin)
2835                  absolutePath=(field[0]=='/');
2836                } else {
2837                  //Windows (non cycgwin)
2838                  absolutePath=(field[0]=='\\');
2839                  // but allow for :
2840                  if (strchr(field.c_str(),':'))
2841                    absolutePath=true;
2842                }
2843                if (absolutePath) {
2844                  fileName = field;
2845                } else if (field[0]=='~') {
2846                  char * environVar = getenv("HOME");
2847                  if (environVar) {
2848                    std::string home(environVar);
2849                    field=field.erase(0,1);
2850                    fileName = home+field;
2851                  } else {
2852                    fileName=field;
2853                  }
2854                } else {
2855                  fileName = directory+field;
2856                }
2857                FILE *fp=fopen(fileName.c_str(),"r");
2858                if (fp) {
2859                  // can open - lets go for it
2860                  fclose(fp);
2861                  canOpen=true;
2862                } else {
2863                  std::cout<<"Unable to open file "<<fileName<<std::endl;
2864                }
2865              }
2866              if (canOpen) {
2867                int status =clpSolver->readMps(fileName.c_str(),
2868                                                   keepImportNames!=0,
2869                                                   allowImportErrors!=0);
2870                if (!status||(status>0&&allowImportErrors)) {
2871                  if (keepImportNames) {
2872                    lengthName = lpSolver->lengthNames();
2873                    rowNames = *(lpSolver->rowNames());
2874                    columnNames = *(lpSolver->columnNames());
2875                  } else {
2876                    lengthName=0;
2877                  }
2878                  goodModel=true;
2879                  //Set integers in clpsolver
2880                  const char * info = lpSolver->integerInformation();
2881                  if (info) {
2882                    int numberColumns = lpSolver->numberColumns();
2883                    int i;
2884                    for (i=0;i<numberColumns;i++) {
2885                      if (info[i]) 
2886                        clpSolver->setInteger(i);
2887                    }
2888                  }
2889                  // sets to all slack (not necessary?)
2890                  lpSolver->createStatus();
2891                  time2 = CoinCpuTime();
2892                  totalTime += time2-time1;
2893                  time1=time2;
2894                  // Go to canned file if just input file
2895                  if (CbcOrClpRead_mode==2&&argc==2) {
2896                    // only if ends .mps
2897                    std::string::size_type loc = fileName.find(".mps") ;
2898                    if (loc != std::string::npos &&
2899                        fileName.length() == loc+3)
2900                    { fileName.replace(loc+1,3,"par") ;
2901                      FILE *fp=fopen(fileName.c_str(),"r");
2902                      if (fp) {
2903                        CbcOrClpReadCommand=fp; // Read from that file
2904                        CbcOrClpRead_mode=-1;
2905                      }
2906                    }
2907                  }
2908                } else {
2909                  // errors
2910                  std::cout<<"There were "<<status<<
2911                    " errors on input"<<std::endl;
2912                }
2913              }
2914            }
2915            break;
2916          case EXPORT:
2917            if (goodModel) {
2918              // get next field
2919              field = CoinReadGetString(argc,argv);
2920              if (field=="$") {
2921                field = parameters[iParam].stringValue();
2922              } else if (field=="EOL") {
2923                parameters[iParam].printString();
2924                break;
2925              } else {
2926                parameters[iParam].setStringValue(field);
2927              }
2928              std::string fileName;
2929              bool canOpen=false;
2930              if (field[0]=='/'||field[0]=='\\') {
2931                fileName = field;
2932              } else if (field[0]=='~') {
2933                char * environVar = getenv("HOME");
2934                if (environVar) {
2935                  std::string home(environVar);
2936                  field=field.erase(0,1);
2937                  fileName = home+field;
2938                } else {
2939                  fileName=field;
2940                }
2941              } else {
2942                fileName = directory+field;
2943              }
2944              FILE *fp=fopen(fileName.c_str(),"w");
2945              if (fp) {
2946                // can open - lets go for it
2947                fclose(fp);
2948                canOpen=true;
2949              } else {
2950                std::cout<<"Unable to open file "<<fileName<<std::endl;
2951              }
2952              if (canOpen) {
2953                // If presolve on then save presolved
2954                bool deleteModel2=false;
2955                ClpSimplex * model2 = lpSolver;
2956#ifdef COIN_HAS_ASL
2957                if (info.numberSos&&doSOS&&usingAmpl) {
2958                  // SOS
2959                  numberSOS = info.numberSos;
2960                  sosStart = info.sosStart;
2961                  sosIndices = info.sosIndices;
2962                  sosReference = info.sosReference;
2963                  preSolve=false;
2964                  clpSolver->setSOSData(numberSOS,info.sosType,sosStart,sosIndices,sosReference);
2965                }
2966#endif
2967                if (preSolve) {
2968                  ClpPresolve pinfo;
2969                  int presolveOptions2 = presolveOptions&~0x40000000;
2970                  if ((presolveOptions2&0xffff)!=0)
2971                    pinfo.setPresolveActions(presolveOptions2);
2972                  if ((printOptions&1)!=0)
2973                    pinfo.statistics();
2974                  double presolveTolerance = 
2975                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
2976                  model2 = 
2977                    pinfo.presolvedModel(*lpSolver,presolveTolerance,
2978                                         true,preSolve);
2979                  if (model2) {
2980                    printf("Saving presolved model on %s\n",
2981                           fileName.c_str());
2982                    deleteModel2=true;
2983                  } else {
2984                    printf("Presolved model looks infeasible - saving original on %s\n",
2985                           fileName.c_str());
2986                    deleteModel2=false;
2987                    model2 = lpSolver;
2988
2989                  }
2990                  model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
2991                  if (deleteModel2)
2992                    delete model2;
2993                } else {
2994                  printf("Saving model on %s\n",
2995                           fileName.c_str());
2996                  if (numberSOS) {
2997                    // Convert names
2998                    int iRow;
2999                    int numberRows=model2->numberRows();
3000                    int iColumn;
3001                    int numberColumns=model2->numberColumns();
3002                   
3003                    char ** rowNames = NULL;
3004                    char ** columnNames = NULL;
3005                    if (model2->lengthNames()) {
3006                      rowNames = new char * [numberRows];
3007                      for (iRow=0;iRow<numberRows;iRow++) {
3008                        rowNames[iRow] = 
3009                          strdup(model2->rowName(iRow).c_str());
3010                      }
3011                     
3012                      columnNames = new char * [numberColumns];
3013                      for (iColumn=0;iColumn<numberColumns;iColumn++) {
3014                        columnNames[iColumn] = 
3015                          strdup(model2->columnName(iColumn).c_str());
3016                      }
3017                    }
3018                    clpSolver->writeMpsNative(fileName.c_str(),(const char **) rowNames,(const char **) columnNames,
3019                                              (outputFormat-1)/2,1+((outputFormat-1)&1));
3020                    if (rowNames) {
3021                      for (iRow=0;iRow<numberRows;iRow++) {
3022                        free(rowNames[iRow]);
3023                      }
3024                      delete [] rowNames;
3025                      for (iColumn=0;iColumn<numberColumns;iColumn++) {
3026                        free(columnNames[iColumn]);
3027                      }
3028                      delete [] columnNames;
3029                    }
3030                  } else {
3031                    model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
3032                  }
3033                }
3034                time2 = CoinCpuTime();
3035                totalTime += time2-time1;
3036                time1=time2;
3037              }
3038            } else {
3039              std::cout<<"** Current model not valid"<<std::endl;
3040            }
3041            break;
3042          case BASISIN:
3043            if (goodModel) {
3044              // get next field
3045              field = CoinReadGetString(argc,argv);
3046              if (field=="$") {
3047                field = parameters[iParam].stringValue();
3048              } else if (field=="EOL") {
3049                parameters[iParam].printString();
3050                break;
3051              } else {
3052                parameters[iParam].setStringValue(field);
3053              }
3054              std::string fileName;
3055              bool canOpen=false;
3056              if (field=="-") {
3057                // stdin
3058                canOpen=true;
3059                fileName = "-";
3060              } else {
3061                if (field[0]=='/'||field[0]=='\\') {
3062                  fileName = field;
3063                } else if (field[0]=='~') {
3064                  char * environVar = getenv("HOME");
3065                  if (environVar) {
3066                    std::string home(environVar);
3067                    field=field.erase(0,1);
3068                    fileName = home+field;
3069                  } else {
3070                    fileName=field;
3071                  }
3072                } else {
3073                  fileName = directory+field;
3074                }
3075                FILE *fp=fopen(fileName.c_str(),"r");
3076                if (fp) {
3077                  // can open - lets go for it
3078                  fclose(fp);
3079                  canOpen=true;
3080                } else {
3081                  std::cout<<"Unable to open file "<<fileName<<std::endl;
3082                }
3083              }
3084              if (canOpen) {
3085                int values = lpSolver->readBasis(fileName.c_str());
3086                if (values==0)
3087                  basisHasValues=-1;
3088                else
3089                  basisHasValues=1;
3090              }
3091            } else {
3092              std::cout<<"** Current model not valid"<<std::endl;
3093            }
3094            break;
3095          case PRIORITYIN:
3096            if (goodModel) {
3097              // get next field
3098              field = CoinReadGetString(argc,argv);
3099              if (field=="$") {
3100                field = parameters[iParam].stringValue();
3101              } else if (field=="EOL") {
3102                parameters[iParam].printString();
3103                break;
3104              } else {
3105                parameters[iParam].setStringValue(field);
3106              }
3107              std::string fileName;
3108              if (field[0]=='/'||field[0]=='\\') {
3109                fileName = field;
3110              } else if (field[0]=='~') {
3111                char * environVar = getenv("HOME");
3112                if (environVar) {
3113                  std::string home(environVar);
3114                  field=field.erase(0,1);
3115                  fileName = home+field;
3116                } else {
3117                  fileName=field;
3118                }
3119              } else {
3120                fileName = directory+field;
3121              }
3122              FILE *fp=fopen(fileName.c_str(),"r");
3123              if (fp) {
3124                // can open - lets go for it
3125                std::string headings[]={"name","number","direction","priority","up","down",
3126                                        "solution","priin"};
3127                int got[]={-1,-1,-1,-1,-1,-1,-1,-1};
3128                int order[8];
3129                assert(sizeof(got)==sizeof(order));
3130                int nAcross=0;
3131                char line[1000];
3132                int numberColumns = lpSolver->numberColumns();
3133                if (!fgets(line,1000,fp)) {
3134                  std::cout<<"Odd file "<<fileName<<std::endl;
3135                } else {
3136                  char * pos = line;
3137                  char * put = line;
3138                  while (*pos>=' '&&*pos!='\n') {
3139                    if (*pos!=' '&&*pos!='\t') {
3140                      *put=tolower(*pos);
3141                      put++;
3142                    }
3143                    pos++;
3144                  }
3145                  *put='\0';
3146                  pos=line;
3147                  int i;
3148                  bool good=true;
3149                  while (pos) {
3150                    char * comma = strchr(pos,',');
3151                    if (comma)
3152                      *comma='\0';
3153                    for (i=0;i<(int) (sizeof(got)/sizeof(int));i++) {
3154                      if (headings[i]==pos) {
3155                        if (got[i]<0) {
3156                          order[nAcross]=i;
3157                          got[i]=nAcross++;
3158                        } else {
3159                          // duplicate
3160                          good=false;
3161                        }
3162                        break;
3163                      }
3164                    }
3165                    if (i==(int) (sizeof(got)/sizeof(int)))
3166                      good=false;
3167                    if (comma) {
3168                      *comma=',';
3169                      pos=comma+1;
3170                    } else {
3171                      break;
3172                    }
3173                  }
3174                  if (got[0]<0&&got[1]<0)
3175                    good=false;
3176                  if (got[0]>=0&&got[1]>=0)
3177                    good=false;
3178                  if (got[0]>=0&&!lpSolver->lengthNames())
3179                    good=false;
3180                  if (good) {
3181                    char ** columnNames = columnNames = new char * [numberColumns];
3182                    pseudoDown= (double *) malloc(numberColumns*sizeof(double));
3183                    pseudoUp = (double *) malloc(numberColumns*sizeof(double));
3184                    branchDirection = (int *) malloc(numberColumns*sizeof(int));
3185                    priorities= (int *) malloc(numberColumns*sizeof(int));
3186                    free(solutionIn);
3187                    solutionIn=NULL;
3188                    free(prioritiesIn);
3189                    prioritiesIn=NULL;
3190                    int iColumn;
3191                    if (got[6]>=0) {
3192                      solutionIn = (double *) malloc(numberColumns*sizeof(double));
3193                      CoinZeroN(solutionIn,numberColumns);
3194                    }
3195                    if (got[7]>=0) {
3196                      prioritiesIn = (int *) malloc(numberColumns*sizeof(int));
3197                      for (iColumn=0;iColumn<numberColumns;iColumn++) 
3198                        prioritiesIn[iColumn]=10000;
3199                    }
3200                    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3201                      columnNames[iColumn] = 
3202                        strdup(lpSolver->columnName(iColumn).c_str());
3203                      pseudoDown[iColumn]=0.0;
3204                      pseudoUp[iColumn]=0.0;
3205                      branchDirection[iColumn]=0;
3206                      priorities[iColumn]=0;
3207                    }
3208                    int nBadPseudo=0;
3209                    int nBadDir=0;
3210                    int nBadPri=0;
3211                    int nBadName=0;
3212                    int nBadLine=0;
3213                    int nLine=0;
3214                    while (fgets(line,1000,fp)) {
3215                      nLine++;
3216                      iColumn = -1;
3217                      double up =0.0;
3218                      double down=0.0;
3219                      int pri=0;
3220                      int dir=0;
3221                      double solValue=COIN_DBL_MAX;
3222                      int priValue=1000000;
3223                      char * pos = line;
3224                      char * put = line;
3225                      while (*pos>=' '&&*pos!='\n') {
3226                        if (*pos!=' '&&*pos!='\t') {
3227                          *put=tolower(*pos);
3228                          put++;
3229                        }
3230                        pos++;
3231                      }
3232                      *put='\0';
3233                      pos=line;
3234                      for (int i=0;i<nAcross;i++) {
3235                        char * comma = strchr(pos,',');
3236                        if (comma) {
3237                          *comma='\0';
3238                        } else if (i<nAcross-1) {
3239                          nBadLine++;
3240                          break;
3241                        }
3242                        switch (order[i]) {
3243                          // name
3244                        case 0:
3245                          for (iColumn=0;iColumn<numberColumns;iColumn++) {
3246                            if (!strcmp(columnNames[iColumn],pos))
3247                              break;
3248                          }
3249                          if (iColumn==numberColumns)
3250                            iColumn=-1;
3251                          break;
3252                          // number
3253                        case 1:
3254                          iColumn = atoi(pos);
3255                          if (iColumn<0||iColumn>=numberColumns)
3256                            iColumn=-1;
3257                          break;
3258                          // direction
3259                        case 2:
3260                          if (*pos=='D')
3261                            dir=-1;
3262                          else if (*pos=='U')
3263                            dir=1;
3264                          else if (*pos=='N')
3265                            dir=0;
3266                          else if (*pos=='1'&&*(pos+1)=='\0')
3267                            dir=1;
3268                          else if (*pos=='0'&&*(pos+1)=='\0')
3269                            dir=0;
3270                          else if (*pos=='1'&&*(pos+1)=='1'&&*(pos+2)=='\0')
3271                            dir=-1;
3272                          else
3273                            dir=-2; // bad
3274                          break;
3275                          // priority
3276                        case 3:
3277                          pri=atoi(pos);
3278                          break;
3279                          // up
3280                        case 4:
3281                          up = atof(pos);
3282                          break;
3283                          // down
3284                        case 5:
3285                          down = atof(pos);
3286                          break;
3287                          // sol value
3288                        case 6:
3289                          solValue = atof(pos);
3290                          break;
3291                          // priority in value
3292                        case 7:
3293                          priValue = atoi(pos);
3294                          break;
3295                        }
3296                        if (comma) {
3297                          *comma=',';
3298                          pos=comma+1;
3299                        }
3300                      }
3301                      if (iColumn>=0) {
3302                        if (down<0.0) {
3303                          nBadPseudo++;
3304                          down=0.0;
3305                        }
3306                        if (up<0.0) {
3307                          nBadPseudo++;
3308                          up=0.0;
3309                        }
3310                        if (!up)
3311                          up=down;
3312                        if (!down)
3313                          down=up;
3314                        if (dir<-1||dir>1) {
3315                          nBadDir++;
3316                          dir=0;
3317                        }
3318                        if (pri<0) {
3319                          nBadPri++;
3320                          pri=0;
3321                        }
3322                        pseudoDown[iColumn]=down;
3323                        pseudoUp[iColumn]=up;
3324                        branchDirection[iColumn]=dir;
3325                        priorities[iColumn]=pri;
3326                        if (solValue!=COIN_DBL_MAX) {
3327                          assert (solutionIn);
3328                          solutionIn[iColumn]=solValue;
3329                        }
3330                        if (priValue!=1000000) {
3331                          assert (prioritiesIn);
3332                          prioritiesIn[iColumn]=priValue;
3333                        }
3334                      } else {
3335                        nBadName++;
3336                      }
3337                    }
3338                    if (!noPrinting) {
3339                      printf("%d fields and %d records",nAcross,nLine);
3340                      if (nBadPseudo)
3341                        printf(" %d bad pseudo costs",nBadPseudo);
3342                      if (nBadDir)
3343                        printf(" %d bad directions",nBadDir);
3344                      if (nBadPri)
3345                        printf(" %d bad priorities",nBadPri);
3346                      if (nBadName)
3347                        printf(" ** %d records did not match on name/sequence",nBadName);
3348                      printf("\n");
3349                    }
3350                    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3351                      free(columnNames[iColumn]);
3352                    }
3353                    delete [] columnNames;
3354                  } else {
3355                    std::cout<<"Duplicate or unknown keyword - or name/number fields wrong"<<line<<std::endl;
3356                  }
3357                }
3358                fclose(fp);
3359              } else {
3360                std::cout<<"Unable to open file "<<fileName<<std::endl;
3361              }
3362            } else {
3363              std::cout<<"** Current model not valid"<<std::endl;
3364            }
3365            break;
3366          case DEBUG:
3367            if (goodModel) {
3368              delete [] debugValues;
3369              debugValues=NULL;
3370              // get next field
3371              field = CoinReadGetString(argc,argv);
3372              if (field=="$") {
3373                field = parameters[iParam].stringValue();
3374              } else if (field=="EOL") {
3375                parameters[iParam].printString();
3376                break;
3377              } else {
3378                parameters[iParam].setStringValue(field);
3379                debugFile=field;
3380                if (debugFile=="create"||
3381                    debugFile=="createAfterPre") {
3382                  printf("Will create a debug file so this run should be a good one\n");
3383                  break;
3384                }
3385              }
3386              std::string fileName;
3387              if (field[0]=='/'||field[0]=='\\') {
3388                fileName = field;
3389              } else if (field[0]=='~') {
3390                char * environVar = getenv("HOME");
3391                if (environVar) {
3392                  std::string home(environVar);
3393                  field=field.erase(0,1);
3394                  fileName = home+field;
3395                } else {
3396                  fileName=field;
3397                }
3398              } else {
3399                fileName = directory+field;
3400              }
3401              FILE *fp=fopen(fileName.c_str(),"rb");
3402              if (fp) {
3403                // can open - lets go for it
3404                int numRows;
3405                double obj;
3406                fread(&numRows,sizeof(int),1,fp);
3407                fread(&numberDebugValues,sizeof(int),1,fp);
3408                fread(&obj,sizeof(double),1,fp);
3409                debugValues = new double[numberDebugValues+numRows];
3410                fread(debugValues,sizeof(double),numRows,fp);
3411                fread(debugValues,sizeof(double),numRows,fp);
3412                fread(debugValues,sizeof(double),numberDebugValues,fp);
3413                printf("%d doubles read into debugValues\n",numberDebugValues);
3414                fclose(fp);
3415              } else {
3416                std::cout<<"Unable to open file "<<fileName<<std::endl;
3417              }
3418            } else {
3419              std::cout<<"** Current model not valid"<<std::endl;
3420            }
3421            break;
3422          case PRINTMASK:
3423            // get next field
3424            {
3425              std::string name = CoinReadGetString(argc,argv);
3426              if (name!="EOL") {
3427                parameters[iParam].setStringValue(name);
3428                printMask = name;
3429              } else {
3430                parameters[iParam].printString();
3431              }
3432            }
3433            break;
3434          case BASISOUT:
3435            if (goodModel) {
3436              // get next field
3437              field = CoinReadGetString(argc,argv);
3438              if (field=="$") {
3439                field = parameters[iParam].stringValue();
3440              } else if (field=="EOL") {
3441                parameters[iParam].printString();
3442                break;
3443              } else {
3444                parameters[iParam].setStringValue(field);
3445              }
3446              std::string fileName;
3447              bool canOpen=false;
3448              if (field[0]=='/'||field[0]=='\\') {
3449                fileName = field;
3450              } else if (field[0]=='~') {
3451                char * environVar = getenv("HOME");
3452                if (environVar) {
3453                  std::string home(environVar);
3454                  field=field.erase(0,1);
3455                  fileName = home+field;
3456                } else {
3457                  fileName=field;
3458                }
3459              } else {
3460                fileName = directory+field;
3461              }
3462              FILE *fp=fopen(fileName.c_str(),"w");
3463              if (fp) {
3464                // can open - lets go for it
3465                fclose(fp);
3466                canOpen=true;
3467              } else {
3468                std::cout<<"Unable to open file "<<fileName<<std::endl;
3469              }
3470              if (canOpen) {
3471                ClpSimplex * model2 = lpSolver;
3472                model2->writeBasis(fileName.c_str(),outputFormat>1,outputFormat-2);
3473                time2 = CoinCpuTime();
3474                totalTime += time2-time1;
3475                time1=time2;
3476              }
3477            } else {
3478              std::cout<<"** Current model not valid"<<std::endl;
3479            }
3480            break;
3481          case SAVE:
3482            {
3483              // get next field
3484              field = CoinReadGetString(argc,argv);
3485              if (field=="$") {
3486                field = parameters[iParam].stringValue();
3487              } else if (field=="EOL") {
3488                parameters[iParam].printString();
3489                break;
3490              } else {
3491                parameters[iParam].setStringValue(field);
3492              }
3493              std::string fileName;
3494              bool canOpen=false;
3495              if (field[0]=='/'||field[0]=='\\') {
3496                fileName = field;
3497              } else if (field[0]=='~') {
3498                char * environVar = getenv("HOME");
3499                if (environVar) {
3500                  std::string home(environVar);
3501                  field=field.erase(0,1);
3502                  fileName = home+field;
3503                } else {
3504                  fileName=field;
3505                }
3506              } else {
3507                fileName = directory+field;
3508              }
3509              FILE *fp=fopen(fileName.c_str(),"wb");
3510              if (fp) {
3511                // can open - lets go for it
3512                fclose(fp);
3513                canOpen=true;
3514              } else {
3515                std::cout<<"Unable to open file "<<fileName<<std::endl;
3516              }
3517              if (canOpen) {
3518                int status;
3519                // If presolve on then save presolved
3520                bool deleteModel2=false;
3521                ClpSimplex * model2 = lpSolver;
3522                if (preSolve) {
3523                  ClpPresolve pinfo;
3524                  double presolveTolerance = 
3525                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
3526                  model2 = 
3527                    pinfo.presolvedModel(*lpSolver,presolveTolerance,
3528                                         false,preSolve);
3529                  if (model2) {
3530                    printf("Saving presolved model on %s\n",
3531                           fileName.c_str());
3532                    deleteModel2=true;
3533                  } else {
3534                    printf("Presolved model looks infeasible - saving original on %s\n",
3535                           fileName.c_str());
3536                    deleteModel2=false;
3537                    model2 = lpSolver;
3538
3539                  }
3540                } else {
3541                  printf("Saving model on %s\n",
3542                           fileName.c_str());
3543                }
3544                status =model2->saveModel(fileName.c_str());
3545                if (deleteModel2)
3546                  delete model2;
3547                if (!status) {
3548                  goodModel=true;
3549                  time2 = CoinCpuTime();
3550                  totalTime += time2-time1;
3551                  time1=time2;
3552                } else {
3553                  // errors
3554                  std::cout<<"There were errors on output"<<std::endl;
3555                }
3556              }
3557            }
3558            break;
3559          case RESTORE:
3560            {
3561              // get next field
3562              field = CoinReadGetString(argc,argv);
3563              if (field=="$") {
3564                field = parameters[iParam].stringValue();
3565              } else if (field=="EOL") {
3566                parameters[iParam].printString();
3567                break;
3568              } else {
3569                parameters[iParam].setStringValue(field);
3570              }
3571              std::string fileName;
3572              bool canOpen=false;
3573              if (field[0]=='/'||field[0]=='\\') {
3574                fileName = field;
3575              } else if (field[0]=='~') {
3576                char * environVar = getenv("HOME");
3577                if (environVar) {
3578                  std::string home(environVar);
3579                  field=field.erase(0,1);
3580                  fileName = home+field;
3581                } else {
3582                  fileName=field;
3583                }
3584              } else {
3585                fileName = directory+field;
3586              }
3587              FILE *fp=fopen(fileName.c_str(),"rb");
3588              if (fp) {
3589                // can open - lets go for it
3590                fclose(fp);
3591                canOpen=true;
3592              } else {
3593                std::cout<<"Unable to open file "<<fileName<<std::endl;
3594              }
3595              if (canOpen) {
3596                int status =lpSolver->restoreModel(fileName.c_str());
3597                if (!status) {
3598                  goodModel=true;
3599                  time2 = CoinCpuTime();
3600                  totalTime += time2-time1;
3601                  time1=time2;
3602                } else {
3603                  // errors
3604                  std::cout<<"There were errors on input"<<std::endl;
3605                }
3606              }
3607            }
3608            break;
3609          case MAXIMIZE:
3610            lpSolver->setOptimizationDirection(-1);
3611            break;
3612          case MINIMIZE:
3613            lpSolver->setOptimizationDirection(1);
3614            break;
3615          case ALLSLACK:
3616            lpSolver->allSlackBasis(true);
3617            break;
3618          case REVERSE:
3619            if (goodModel) {
3620              int iColumn;
3621              int numberColumns=lpSolver->numberColumns();
3622              double * dualColumnSolution = 
3623                lpSolver->dualColumnSolution();
3624              ClpObjective * obj = lpSolver->objectiveAsObject();
3625              assert(dynamic_cast<ClpLinearObjective *> (obj));
3626              double offset;
3627              double * objective = obj->gradient(NULL,NULL,offset,true);
3628              for (iColumn=0;iColumn<numberColumns;iColumn++) {
3629                dualColumnSolution[iColumn] = dualColumnSolution[iColumn];
3630                objective[iColumn] = -objective[iColumn];
3631              }
3632              int iRow;
3633              int numberRows=lpSolver->numberRows();
3634              double * dualRowSolution = 
3635                lpSolver->dualRowSolution();
3636              for (iRow=0;iRow<numberRows;iRow++) 
3637                dualRowSolution[iRow] = dualRowSolution[iRow];
3638            }
3639            break;
3640          case DIRECTORY:
3641            {
3642              std::string name = CoinReadGetString(argc,argv);
3643              if (name!="EOL") {
3644                int length=name.length();
3645                if (name[length-1]=='/'||name[length-1]=='\\')
3646                  directory=name;
3647                else
3648                  directory = name+"/";
3649                parameters[iParam].setStringValue(directory);
3650              } else {
3651                parameters[iParam].printString();
3652              }
3653            }
3654            break;
3655          case STDIN:
3656            CbcOrClpRead_mode=-1;
3657            break;
3658          case NETLIB_DUAL:
3659          case NETLIB_EITHER:
3660          case NETLIB_BARRIER:
3661          case NETLIB_PRIMAL:
3662          case NETLIB_TUNE:
3663            {
3664              // create fields for unitTest
3665              const char * fields[4];
3666              int nFields=2;
3667              fields[0]="fake main from unitTest";
3668              fields[1]="-netlib";
3669              if (directory!=defaultDirectory) {
3670                fields[2]="-netlibDir";
3671                fields[3]=directory.c_str();
3672                nFields=4;
3673              }
3674              int algorithm;
3675              if (type==NETLIB_DUAL) {
3676                std::cerr<<"Doing netlib with dual algorithm"<<std::endl;
3677                algorithm =0;
3678              } else if (type==NETLIB_BARRIER) {
3679                std::cerr<<"Doing netlib with barrier algorithm"<<std::endl;
3680                algorithm =2;
3681              } else if (type==NETLIB_EITHER) {
3682                std::cerr<<"Doing netlib with dual or primal algorithm"<<std::endl;
3683                algorithm =3;
3684              } else if (type==NETLIB_TUNE) {
3685                std::cerr<<"Doing netlib with best algorithm!"<<std::endl;
3686                algorithm =5;
3687                // uncomment next to get active tuning
3688                // algorithm=6;
3689              } else {
3690                std::cerr<<"Doing netlib with primal agorithm"<<std::endl;
3691                algorithm=1;
3692              }
3693              int specialOptions = lpSolver->specialOptions();
3694              lpSolver->setSpecialOptions(0);
3695              mainTest(nFields,fields,algorithm,*lpSolver,
3696                       (preSolve!=0),specialOptions,doVector!=0);
3697            }
3698            break;
3699          case UNITTEST:
3700            {
3701              // create fields for unitTest
3702              const char * fields[3];
3703              int nFields=1;
3704              fields[0]="fake main from unitTest";
3705              if (directory!=defaultDirectory) {
3706                fields[1]="-mpsDir";
3707                fields[2]=directory.c_str();
3708                nFields=3;
3709              }
3710              mainTest(nFields,fields,false,*lpSolver,(preSolve!=0),
3711                       false,doVector!=0);
3712            }
3713            break;
3714          case FAKEBOUND:
3715            if (goodModel) {
3716              // get bound
3717              double value = CoinReadGetDoubleField(argc,argv,&valid);
3718              if (!valid) {
3719                std::cout<<"Setting "<<parameters[iParam].name()<<
3720                  " to DEBUG "<<value<<std::endl;
3721                int iRow;
3722                int numberRows=lpSolver->numberRows();
3723                double * rowLower = lpSolver->rowLower();
3724                double * rowUpper = lpSolver->rowUpper();
3725                for (iRow=0;iRow<numberRows;iRow++) {
3726                  // leave free ones for now
3727                  if (rowLower[iRow]>-1.0e20||rowUpper[iRow]<1.0e20) {
3728                    rowLower[iRow]=CoinMax(rowLower[iRow],-value);
3729                    rowUpper[iRow]=CoinMin(rowUpper[iRow],value);
3730                  }
3731                }
3732                int iColumn;
3733                int numberColumns=lpSolver->numberColumns();
3734                double * columnLower = lpSolver->columnLower();
3735                double * columnUpper = lpSolver->columnUpper();
3736                for (iColumn=0;iColumn<numberColumns;iColumn++) {
3737                  // leave free ones for now
3738                  if (columnLower[iColumn]>-1.0e20||
3739                      columnUpper[iColumn]<1.0e20) {
3740                    columnLower[iColumn]=CoinMax(columnLower[iColumn],-value);
3741                    columnUpper[iColumn]=CoinMin(columnUpper[iColumn],value);
3742                  }
3743                }
3744              } else if (valid==1) {
3745                abort();
3746              } else {
3747                std::cout<<"enter value for "<<parameters[iParam].name()<<
3748                  std::endl;
3749              }
3750            }
3751            break;
3752          case REALLY_SCALE:
3753            if (goodModel) {
3754              ClpSimplex newModel(*lpSolver,
3755                                  lpSolver->scalingFlag());
3756              printf("model really really scaled\n");
3757              *lpSolver=newModel;
3758            }
3759            break;
3760          case USERCLP:
3761            // Replace the sample code by whatever you want
3762            if (goodModel) {
3763              printf("Dummy user clp code - model has %d rows and %d columns\n",
3764                     lpSolver->numberRows(),lpSolver->numberColumns());
3765            }
3766            break;
3767          case USERCBC:
3768            // Replace the sample code by whatever you want
3769            if (goodModel) {
3770#ifndef USER_HAS_FAKE_MAIN
3771              printf("Dummy user cbc code - model has %d rows and %d columns\n",
3772                     model.getNumRows(),model.getNumCols());
3773              // Reduce printout
3774              model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
3775              // Do complete search
3776              model.branchAndBound();
3777#ifdef COIN_HAS_ASL
3778              double objectiveValue=model.getMinimizationObjValue();
3779              int iStat = model.status();
3780              int iStat2 = model.secondaryStatus();
3781#endif
3782#else
3783              // Way of using an existing piece of code
3784              OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
3785              ClpSimplex * lpSolver = clpSolver->getModelPtr();
3786              // set time from integer model
3787              double timeToGo = model.getMaximumSeconds();
3788              lpSolver->setMaximumSeconds(timeToGo);
3789              fakeMain(*lpSolver,*clpSolver,model);
3790#ifdef COIN_HAS_ASL
3791              // My actual usage has objective only in clpSolver
3792              double objectiveValue=clpSolver->getObjValue();
3793              int iStat = lpSolver->status();
3794              int iStat2 = lpSolver->secondaryStatus();
3795#endif
3796#endif
3797              // make sure solution back in correct place
3798              clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
3799              lpSolver = clpSolver->getModelPtr();
3800#ifdef COIN_HAS_ASL
3801              if (usingAmpl) {
3802                int n = clpSolver->getNumCols();
3803                double value = objectiveValue*lpSolver->getObjSense();
3804                char buf[300];
3805                int pos=0;
3806                std::string minor[]={"","","gap","nodes","time","","solutions","user ctrl-c"};
3807                if (iStat==0) {
3808                  if (objectiveValue<1.0e40) {
3809                    pos += sprintf(buf+pos,"optimal," );
3810                  } else {
3811                    // infeasible
3812                    iStat=1;
3813                    pos += sprintf(buf+pos,"infeasible,");
3814                  }
3815                } else if (iStat==1) {
3816                  if (iStat2!=6)
3817                    iStat=3;
3818                  else
3819                    iStat=4;
3820                  pos += sprintf(buf+pos,"stopped on %s,",minor[iStat2].c_str());
3821                } else if (iStat==2) {
3822                  iStat = 7;
3823                  pos += sprintf(buf+pos,"stopped on difficulties,");
3824                } else if (iStat==5) {
3825                  iStat = 3;
3826                  pos += sprintf(buf+pos,"stopped on ctrl-c,");
3827                } else {
3828                  pos += sprintf(buf+pos,"status unknown,");
3829                  iStat=6;
3830                }
3831                info.problemStatus=iStat;
3832                info.objValue = value;
3833                if (objectiveValue<1.0e40) 
3834                  pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
3835                                 value);
3836                sprintf(buf+pos,"\n%d nodes, %d iterations",
3837                        model.getNodeCount(),
3838                        model.getIterationCount());
3839                if (objectiveValue<1.0e50) {
3840                  free(info.primalSolution);
3841                  info.primalSolution = (double *) malloc(n*sizeof(double));
3842                  CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution);
3843                  int numberRows = lpSolver->numberRows();
3844                  free(info.dualSolution);
3845                  info.dualSolution = (double *) malloc(numberRows*sizeof(double));
3846                  CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution);
3847                } else {
3848                  info.primalSolution=NULL;
3849                  info.dualSolution=NULL;
3850                }
3851                // put buffer into info
3852                strcpy(info.buffer,buf);
3853              }
3854#endif
3855            }
3856            break;
3857          case HELP:
3858            std::cout<<"Coin Solver version "<<CBCVERSION
3859                     <<", build "<<__DATE__<<std::endl;
3860            std::cout<<"Non default values:-"<<std::endl;
3861            std::cout<<"Perturbation "<<lpSolver->perturbation()<<" (default 100)"
3862                     <<std::endl;
3863            CoinReadPrintit(
3864                    "Presolve being done with 5 passes\n\
3865Dual steepest edge steep/partial on matrix shape and factorization density\n\
3866Clpnnnn taken out of messages\n\
3867If Factorization frequency default then done on size of matrix\n\n\
3868(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
3869You can switch to interactive mode at any time so\n\
3870clp watson.mps -scaling off -primalsimplex\nis the same as\n\
3871clp watson.mps -\nscaling off\nprimalsimplex"
3872                    );
3873            break;
3874          case SOLUTION:
3875            if (goodModel) {
3876              // get next field
3877              field = CoinReadGetString(argc,argv);
3878              if (field=="$") {
3879                field = parameters[iParam].stringValue();
3880              } else if (field=="EOL") {
3881                parameters[iParam].printString();
3882                break;
3883              } else {
3884                parameters[iParam].setStringValue(field);
3885              }
3886              std::string fileName;
3887              FILE *fp=NULL;
3888              if (field=="-"||field=="EOL"||field=="stdout") {
3889                // stdout
3890                fp=stdout;
3891              } else if (field=="stderr") {
3892                // stderr
3893                fp=stderr;
3894              } else {
3895                if (field[0]=='/'||field[0]=='\\') {
3896                  fileName = field;
3897                } else if (field[0]=='~') {
3898                  char * environVar = getenv("HOME");
3899                  if (environVar) {
3900                    std::string home(environVar);
3901                    field=field.erase(0,1);
3902                    fileName = home+field;
3903                  } else {
3904                    fileName=field;
3905                  }
3906                } else {
3907                  fileName = directory+field;
3908                }
3909                fp=fopen(fileName.c_str(),"w");
3910              }
3911              if (fp) {
3912                // make fancy later on
3913                int iRow;
3914                int numberRows=lpSolver->numberRows();
3915                double * dualRowSolution = lpSolver->dualRowSolution();
3916                double * primalRowSolution = 
3917                  lpSolver->primalRowSolution();
3918                double * rowLower = lpSolver->rowLower();
3919                double * rowUpper = lpSolver->rowUpper();
3920                double primalTolerance = lpSolver->primalTolerance();
3921                char format[6];
3922                sprintf(format,"%%-%ds",CoinMax(lengthName,8));
3923                bool doMask = (printMask!=""&&lengthName);
3924                int * maskStarts=NULL;
3925                int maxMasks=0;
3926                char ** masks =NULL;
3927                if (doMask) {
3928                  int nAst =0;
3929                  const char * pMask2 = printMask.c_str();
3930                  char pMask[100];
3931                  int iChar;
3932                  int lengthMask = strlen(pMask2);
3933                  assert (lengthMask<100);
3934                  if (*pMask2=='"') {
3935                    if (pMask2[lengthMask-1]!='"') {
3936                      printf("mismatched \" in mask %s\n",pMask2);
3937                      break;
3938                    } else {
3939                      strcpy(pMask,pMask2+1);
3940                      *strchr(pMask,'"')='\0';
3941                    }
3942                  } else if (*pMask2=='\'') {
3943                    if (pMask2[lengthMask-1]!='\'') {
3944                      printf("mismatched ' in mask %s\n",pMask2);
3945                      break;
3946                    } else {
3947                      strcpy(pMask,pMask2+1);
3948                      *strchr(pMask,'\'')='\0';
3949                    }
3950                  } else {
3951                    strcpy(pMask,pMask2);
3952                  }
3953                  if (lengthMask>lengthName) {
3954                    printf("mask %s too long - skipping\n",pMask);
3955                    break;
3956                  }
3957                  maxMasks = 1;
3958                  for (iChar=0;iChar<lengthMask;iChar++) {
3959                    if (pMask[iChar]=='*') {
3960                      nAst++;
3961                      maxMasks *= (lengthName+1);
3962                    }
3963                  }
3964                  int nEntries = 1;
3965                  maskStarts = new int[lengthName+2];
3966                  masks = new char * [maxMasks];
3967                  char ** newMasks = new char * [maxMasks];
3968                  int i;
3969                  for (i=0;i<maxMasks;i++) {
3970                    masks[i] = new char[lengthName+1];
3971                    newMasks[i] = new char[lengthName+1];
3972                  }
3973                  strcpy(masks[0],pMask);
3974                  for (int iAst=0;iAst<nAst;iAst++) {
3975                    int nOldEntries = nEntries;
3976                    nEntries=0;
3977                    for (int iEntry = 0;iEntry<nOldEntries;iEntry++) {
3978                      char * oldMask = masks[iEntry];
3979                      char * ast = strchr(oldMask,'*');
3980                      assert (ast);
3981                      int length = strlen(oldMask)-1;
3982                      int nBefore = ast-oldMask;
3983                      int nAfter = length-nBefore;
3984                      // and add null
3985                      nAfter++;
3986                      for (int i=0;i<=lengthName-length;i++) {
3987                        char * maskOut = newMasks[nEntries];
3988                        memcpy(maskOut,oldMask,nBefore);
3989                        for (int k=0;k<i;k++) 
3990                          maskOut[k+nBefore]='?';
3991                        memcpy(maskOut+nBefore+i,ast+1,nAfter);
3992                        nEntries++;
3993                        assert (nEntries<=maxMasks);
3994                      }
3995                    }
3996                    char ** temp = masks;
3997                    masks = newMasks;
3998                    newMasks = temp;
3999                  }
4000                  // Now extend and sort
4001                  int * sort = new int[nEntries];
4002                  for (i=0;i<nEntries;i++) {
4003                    char * maskThis = masks[i];
4004                    int length = strlen(maskThis);
4005                    while (maskThis[length-1]==' ')
4006                      length--;
4007                    maskThis[length]='\0';
4008                    sort[i]=length;
4009                  }
4010                  CoinSort_2(sort,sort+nEntries,masks);
4011                  int lastLength=-1;
4012                  for (i=0;i<nEntries;i++) {
4013                    int length = sort[i];
4014                    while (length>lastLength) 
4015                      maskStarts[++lastLength] = i;
4016                  }
4017                  maskStarts[++lastLength]=nEntries;
4018                  delete [] sort;
4019                  for (i=0;i<maxMasks;i++)
4020                    delete [] newMasks[i];
4021                  delete [] newMasks;
4022                }
4023                if (printMode>2) {
4024                  for (iRow=0;iRow<numberRows;iRow++) {
4025                    int type=printMode-3;
4026                    if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
4027                        primalRowSolution[iRow]<rowLower[iRow]-primalTolerance) {
4028                      fprintf(fp,"** ");
4029                      type=2;
4030                    } else if (fabs(primalRowSolution[iRow])>1.0e-8) {
4031                      type=1;
4032                    } else if (numberRows<50) {
4033                      type=3;
4034                    }
4035                    if (doMask&&!maskMatches(maskStarts,masks,rowNames[iRow]))
4036                      type =0;
4037                    if (type) {
4038                      fprintf(fp,"%7d ",iRow);
4039                      if (lengthName)
4040                        fprintf(fp,format,rowNames[iRow].c_str());
4041                      fprintf(fp,"%15.8g        %15.8g\n",primalRowSolution[iRow],
4042                              dualRowSolution[iRow]);
4043                    }
4044                  }
4045                }
4046                int iColumn;
4047                int numberColumns=lpSolver->numberColumns();
4048                double * dualColumnSolution = 
4049                  lpSolver->dualColumnSolution();
4050                double * primalColumnSolution = 
4051                  lpSolver->primalColumnSolution();
4052                double * columnLower = lpSolver->columnLower();
4053                double * columnUpper = lpSolver->columnUpper();
4054                if (printMode!=2) {
4055                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4056                    int type=(printMode>3) ? 1 :0;
4057                    if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
4058                        primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance) {
4059                      fprintf(fp,"** ");
4060                      type=2;
4061                    } else if (fabs(primalColumnSolution[iColumn])>1.0e-8) {
4062                      type=1;
4063                    } else if (numberColumns<50) {
4064                      type=3;
4065                    }
4066                    // see if integer
4067                    if ((!lpSolver->isInteger(iColumn)||fabs(primalColumnSolution[iColumn])<1.0e-8)
4068                         &&printMode==1)
4069                      type=0;
4070                    if (doMask&&!maskMatches(maskStarts,masks,
4071                                             columnNames[iColumn]))
4072                      type =0;
4073                    if (type) {
4074                      fprintf(fp,"%7d ",iColumn);
4075                      if (lengthName)
4076                        fprintf(fp,format,columnNames[iColumn].c_str());
4077                      fprintf(fp,"%15.8g        %15.8g\n",
4078                              primalColumnSolution[iColumn],
4079                              dualColumnSolution[iColumn]);
4080                    }
4081                  }
4082                } else {
4083                  // special format suitable for OsiRowCutDebugger
4084                  int n=0;
4085                  bool comma=false;
4086                  bool newLine=false;
4087                  fprintf(fp,"\tint intIndicesV[]={\n");
4088                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4089                    if(primalColumnSolution[iColumn]>0.5&&model.solver()->isInteger(iColumn)) {
4090                      if (comma)
4091                        fprintf(fp,",");
4092                      if (newLine)
4093                        fprintf(fp,"\n");
4094                      fprintf(fp,"%d ",iColumn);
4095                      comma=true;
4096                      newLine=false;
4097                      n++;
4098                      if (n==10) {
4099                        n=0;
4100                        newLine=true;
4101                      }
4102                    }
4103                  }
4104                  fprintf(fp,"};\n");
4105                  n=0;
4106                  comma=false;
4107                  newLine=false;
4108                  fprintf(fp,"\tdouble intSolnV[]={\n");
4109                  for ( iColumn=0;iColumn<numberColumns;iColumn++) {
4110                    if(primalColumnSolution[iColumn]>0.5&&model.solver()->isInteger(iColumn)) {
4111                      if (comma)
4112                        fprintf(fp,",");
4113                      if (newLine)
4114                        fprintf(fp,"\n");
4115                      int value = (int) (primalColumnSolution[iColumn]+0.5);
4116                      fprintf(fp,"%d. ",value);
4117                      comma=true;
4118                      newLine=false;
4119                      n++;
4120                      if (n==10) {
4121                        n=0;
4122                        newLine=true;
4123                      }
4124                    }
4125                  }
4126                  fprintf(fp,"};\n");
4127                }
4128                if (fp!=stdout)
4129                  fclose(fp);
4130                if (masks) {
4131                  delete [] maskStarts;
4132                  for (int i=0;i<maxMasks;i++)
4133                    delete [] masks[i];
4134                  delete [] masks;
4135                }
4136              } else {
4137                std::cout<<"Unable to open file "<<fileName<<std::endl;
4138              }
4139            } else {
4140              std::cout<<"** Current model not valid"<<std::endl;
4141             
4142            }
4143            break;
4144          case SAVESOL:
4145            if (goodModel) {
4146              // get next field
4147              field = CoinReadGetString(argc,argv);
4148              if (field=="$") {
4149                field = parameters[iParam].stringValue();
4150              } else if (field=="EOL") {
4151                parameters[iParam].printString();
4152                break;
4153              } else {
4154                parameters[iParam].setStringValue(field);
4155              }
4156              std::string fileName;
4157              if (field[0]=='/'||field[0]=='\\') {
4158                fileName = field;
4159              } else if (field[0]=='~') {
4160                char * environVar = getenv("HOME");
4161                if (environVar) {
4162                  std::string home(environVar);
4163                  field=field.erase(0,1);
4164                  fileName = home+field;
4165                } else {
4166                  fileName=field;
4167                }
4168              } else {
4169                fileName = directory+field;
4170              }
4171              saveSolution(lpSolver,fileName);
4172            } else {
4173              std::cout<<"** Current model not valid"<<std::endl;
4174             
4175            }
4176            break;
4177          case DUMMY:
4178            break;
4179          default:
4180            abort();
4181          }
4182        } 
4183      } else if (!numberMatches) {
4184        std::cout<<"No match for "<<field<<" - ? for list of commands"
4185                 <<std::endl;
4186      } else if (numberMatches==1) {
4187        if (!numberQuery) {
4188          std::cout<<"Short match for "<<field<<" - completion: ";
4189          std::cout<<parameters[firstMatch].matchName()<<std::endl;
4190        } else if (numberQuery) {
4191          std::cout<<parameters[firstMatch].matchName()<<" : ";
4192          std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
4193          if (numberQuery>=2) 
4194            parameters[firstMatch].printLongHelp();
4195        }
4196      } else {
4197        if (!numberQuery) 
4198          std::cout<<"Multiple matches for "<<field<<" - possible completions:"
4199                   <<std::endl;
4200        else
4201          std::cout<<"Completions of "<<field<<":"<<std::endl;
4202        for ( iParam=0; iParam<numberParameters; iParam++ ) {
4203          int match = parameters[iParam].matches(field);
4204          if (match&&parameters[iParam].displayThis()) {
4205            std::cout<<parameters[iParam].matchName();
4206            if (numberQuery>=2) 
4207              std::cout<<" : "<<parameters[iParam].shortHelp();
4208            std::cout<<std::endl;
4209          }
4210        }
4211      }
4212    }
4213  }
4214  // By now all memory should be freed
4215#ifdef DMALLOC
4216  dmalloc_log_unfreed();
4217  dmalloc_shutdown();
4218#endif
4219  return 0;
4220}   
4221static void breakdown(const char * name, int numberLook, const double * region)
4222{
4223  double range[] = {
4224    -COIN_DBL_MAX,
4225    -1.0e15,-1.0e11,-1.0e8,-1.0e5,-1.0e4,-1.0e3,-1.0e2,-1.0e1,
4226    -1.0,
4227    -1.0e-1,-1.0e-2,-1.0e-3,-1.0e-4,-1.0e-5,-1.0e-8,-1.0e-11,-1.0e-15,
4228    0.0,
4229    1.0e-15,1.0e-11,1.0e-8,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,
4230    1.0,
4231    1.0e1,1.0e2,1.0e3,1.0e4,1.0e5,1.0e8,1.0e11,1.0e15,
4232    COIN_DBL_MAX};
4233  int nRanges = (int) (sizeof(range)/sizeof(double));
4234  int * number = new int[nRanges];
4235  memset(number,0,nRanges*sizeof(int));
4236  int * numberExact = new int[nRanges];
4237  memset(numberExact,0,nRanges*sizeof(int));
4238  int i;
4239  for ( i=0;i<numberLook;i++) {
4240    double value = region[i];
4241    for (int j=0;j<nRanges;j++) {
4242      if (value==range[j]) {
4243        numberExact[j]++;
4244        break;
4245      } else if (value<range[j]) {
4246        number[j]++;
4247        break;
4248      }
4249    }
4250  }
4251  printf("\n%s has %d entries\n",name,numberLook);
4252  for (i=0;i<nRanges;i++) {
4253    if (number[i]) 
4254      printf("%d between %g and %g",number[i],range[i-1],range[i]);
4255    if (numberExact[i]) {
4256      if (number[i])
4257        printf(", ");
4258      printf("%d exactly at %g",numberExact[i],range[i]);
4259    }
4260    if (number[i]+numberExact[i])
4261      printf("\n");
4262  }
4263  delete [] number;
4264  delete [] numberExact;
4265}
4266static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
4267{
4268  int numberColumns = originalModel->numberColumns();
4269  const char * integerInformation  = originalModel->integerInformation(); 
4270  const double * columnLower = originalModel->columnLower();
4271  const double * columnUpper = originalModel->columnUpper();
4272  int numberIntegers=0;
4273  int numberBinary=0;
4274  int iRow,iColumn;
4275  if (integerInformation) {
4276    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4277      if (integerInformation[iColumn]) {
4278        if (columnUpper[iColumn]>columnLower[iColumn]) {
4279          numberIntegers++;
4280          if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==1) 
4281            numberBinary++;
4282        }
4283      }
4284    }
4285  }
4286  numberColumns = model->numberColumns();
4287  int numberRows = model->numberRows();
4288  columnLower = model->columnLower();
4289  columnUpper = model->columnUpper();
4290  const double * rowLower = model->rowLower();
4291  const double * rowUpper = model->rowUpper();
4292  const double * objective = model->objective();
4293  CoinPackedMatrix * matrix = model->matrix();
4294  CoinBigIndex numberElements = matrix->getNumElements();
4295  const int * columnLength = matrix->getVectorLengths();
4296  //const CoinBigIndex * columnStart = matrix->getVectorStarts();
4297  const double * elementByColumn = matrix->getElements();
4298  int * number = new int[numberRows+1];
4299  memset(number,0,(numberRows+1)*sizeof(int));
4300  int numberObjSingletons=0;
4301  /* cType
4302     0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
4303     8 0/1
4304  */ 
4305  int cType[9];
4306  std::string cName[]={"0.0->inf,","0.0->up,","lo->inf,","lo->up,","free,","fixed,","-inf->0.0,",
4307                       "-inf->up,","0.0->1.0"};
4308  int nObjective=0;
4309  memset(cType,0,sizeof(cType));
4310  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4311    int length=columnLength[iColumn];
4312    if (length==1&&objective[iColumn])
4313      numberObjSingletons++;
4314    number[length]++;
4315    if (objective[iColumn])
4316      nObjective++;
4317    if (columnLower[iColumn]>-1.0e20) {
4318      if (columnLower[iColumn]==0.0) {
4319        if (columnUpper[iColumn]>1.0e20)
4320          cType[0]++;
4321        else if (columnUpper[iColumn]==1.0)
4322          cType[8]++;
4323        else if (columnUpper[iColumn]==0.0)
4324          cType[5]++;
4325        else
4326          cType[1]++;
4327      } else {
4328        if (columnUpper[iColumn]>1.0e20) 
4329          cType[2]++;
4330        else if (columnUpper[iColumn]==columnLower[iColumn])
4331          cType[5]++;
4332        else
4333          cType[3]++;
4334      }
4335    } else {
4336      if (columnUpper[iColumn]>1.0e20) 
4337        cType[4]++;
4338      else if (columnUpper[iColumn]==0.0) 
4339        cType[6]++;
4340      else
4341        cType[7]++;
4342    }
4343  }
4344  /* rType
4345     0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
4346     7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
4347  */ 
4348  int rType[13];
4349  std::string rName[]={"E 0.0,","E 1.0,","E -1.0,","E other,","G 0.0,","G 1.0,","G other,",
4350                       "L 0.0,","L 1.0,","L other,","Range 0.0->1.0,","Range other,","Free"};
4351  memset(rType,0,sizeof(rType));
4352  for (iRow=0;iRow<numberRows;iRow++) {
4353    if (rowLower[iRow]>-1.0e20) {
4354      if (rowLower[iRow]==0.0) {
4355        if (rowUpper[iRow]>1.0e20)
4356          rType[4]++;
4357        else if (rowUpper[iRow]==1.0)
4358          rType[10]++;
4359        else if (rowUpper[iRow]==0.0)
4360          rType[0]++;
4361        else
4362          rType[11]++;
4363      } else if (rowLower[iRow]==1.0) {
4364        if (rowUpper[iRow]>1.0e20) 
4365          rType[5]++;
4366        else if (rowUpper[iRow]==rowLower[iRow])
4367          rType[1]++;
4368        else
4369          rType[11]++;
4370      } else if (rowLower[iRow]==-1.0) {
4371        if (rowUpper[iRow]>1.0e20) 
4372          rType[6]++;
4373        else if (rowUpper[iRow]==rowLower[iRow])
4374          rType[2]++;
4375        else
4376          rType[11]++;
4377      } else {
4378        if (rowUpper[iRow]>1.0e20) 
4379          rType[6]++;
4380        else if (rowUpper[iRow]==rowLower[iRow])
4381          rType[3]++;
4382        else
4383          rType[11]++;
4384      }
4385    } else {
4386      if (rowUpper[iRow]>1.0e20) 
4387        rType[12]++;
4388      else if (rowUpper[iRow]==0.0) 
4389        rType[7]++;
4390      else if (rowUpper[iRow]==1.0) 
4391        rType[8]++;
4392      else
4393        rType[9]++;
4394    }
4395  }
4396  // Basic statistics
4397  printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
4398         numberRows,numberColumns,nObjective,numberElements);
4399  if (number[0]+number[1]) {
4400    printf("There are ");
4401    if (numberObjSingletons)
4402      printf("%d singletons with objective ",numberObjSingletons);
4403    int numberNoObj = number[1]-numberObjSingletons;
4404    if (numberNoObj)
4405      printf("%d singletons with no objective ",numberNoObj);
4406    if (number[0])
4407      printf("** %d columns have no entries",number[0]);
4408    printf("\n");
4409  }
4410  printf("Column breakdown:\n");
4411  int k;
4412  for (k=0;k<(int) (sizeof(cType)/sizeof(int));k++) {
4413    printf("%d of type %s ",cType[k],cName[k].c_str());
4414    if (((k+1)%3)==0)
4415      printf("\n");
4416  }
4417  if ((k%3)!=0)
4418    printf("\n");
4419  printf("Row breakdown:\n");
4420  for (k=0;k<(int) (sizeof(rType)/sizeof(int));k++) {
4421    printf("%d of type %s ",rType[k],rName[k].c_str());
4422    if (((k+1)%3)==0)
4423      printf("\n");
4424  }
4425  if ((k%3)!=0)
4426    printf("\n");
4427  if (model->logLevel()<2)
4428    return ;
4429  int kMax = model->logLevel()>3 ? 1000000 : 10;
4430  k=0;
4431  for (iRow=1;iRow<=numberRows;iRow++) {
4432    if (number[iRow]) {
4433      k++;
4434      printf("%d columns have %d entries\n",number[iRow],iRow);
4435      if (k==kMax)
4436        break;
4437    }
4438  }
4439  if (k<numberRows) {
4440    int kk=k;
4441    k=0;
4442    for (iRow=numberRows;iRow>=1;iRow--) {
4443      if (number[iRow]) {
4444        k++;
4445        if (k==kMax)
4446          break;
4447      }
4448    }
4449    if (k>kk) {
4450      printf("\n    .........\n\n");
4451      iRow=k;
4452      k=0;
4453      for (;iRow<numberRows;iRow++) {
4454        if (number[iRow]) {
4455          k++;
4456          printf("%d columns have %d entries\n",number[iRow],iRow);
4457          if (k==kMax)
4458            break;
4459        }
4460      }
4461    }
4462  }
4463  delete [] number;
4464  printf("\n\n");
4465  // get row copy
4466  CoinPackedMatrix rowCopy = *matrix;
4467  rowCopy.reverseOrdering();
4468  //const int * column = rowCopy.getIndices();
4469  const int * rowLength = rowCopy.getVectorLengths();
4470  //const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
4471  //const double * element = rowCopy.getElements();
4472  number = new int[numberColumns+1];
4473  memset(number,0,(numberColumns+1)*sizeof(int));
4474  for (iRow=0;iRow<numberRows;iRow++) {
4475    int length=rowLength[iRow];
4476    number[length]++;
4477  }
4478  if (number[0])
4479    printf("** %d rows have no entries\n",number[0]);
4480  k=0;
4481  for (iColumn=1;iColumn<=numberColumns;iColumn++) {
4482    if (number[iColumn]) {
4483      k++;
4484      printf("%d rows have %d entries\n",number[iColumn],iColumn);
4485      if (k==kMax)
4486        break;
4487    }
4488  }
4489  if (k<numberColumns) {
4490    int kk=k;
4491    k=0;
4492    for (iColumn=numberColumns;iColumn>=1;iColumn--) {
4493      if (number[iColumn]) {
4494        k++;
4495        if (k==kMax)
4496          break;
4497      }
4498    }
4499    if (k>kk) {
4500      printf("\n    .........\n\n");
4501      iColumn=k;
4502      k=0;
4503      for (;iColumn<numberColumns;iColumn++) {
4504        if (number[iColumn]) {
4505          k++;
4506          printf("%d rows have %d entries\n",number[iColumn],iColumn);
4507          if (k==kMax)
4508            break;
4509        }
4510      }
4511    }
4512  }
4513  delete [] number;
4514  // Now do breakdown of ranges
4515  breakdown("Elements",numberElements,elementByColumn);
4516  breakdown("RowLower",numberRows,rowLower);
4517  breakdown("RowUpper",numberRows,rowUpper);
4518  breakdown("ColumnLower",numberColumns,columnLower);
4519  breakdown("ColumnUpper",numberColumns,columnUpper);
4520  breakdown("Objective",numberColumns,objective);
4521}
4522static bool maskMatches(const int * starts, char ** masks,
4523                        std::string & check)
4524{
4525  // back to char as I am old fashioned
4526  const char * checkC = check.c_str();
4527  int length = strlen(checkC);
4528  while (checkC[length-1]==' ')
4529    length--;
4530  for (int i=starts[length];i<starts[length+1];i++) {
4531    char * thisMask = masks[i];
4532    int k;
4533    for ( k=0;k<length;k++) {
4534      if (thisMask[k]!='?'&&thisMask[k]!=checkC[k]) 
4535        break;
4536    }
4537    if (k==length)
4538      return true;
4539  }
4540  return false;
4541}
4542static void clean(char * temp)
4543{
4544  char * put = temp;
4545  while (*put>=' ')
4546    put++;
4547  *put='\0';
4548}
4549static void generateCode(CbcModel * model, const char * fileName,int type,int preProcess)
4550{
4551  // options on code generation
4552  bool sizecode = (type&4)!=0;
4553  type &= 3;
4554  FILE * fp = fopen(fileName,"r");
4555  assert (fp);
4556  int numberLines=0;
4557#define MAXLINES 5000
4558#define MAXONELINE 200
4559  char line[MAXLINES][MAXONELINE];
4560  strcpy(line[numberLines++],"0#if defined(_MSC_VER)");
4561  strcpy(line[numberLines++],"0// Turn off compiler warning about long names");
4562  strcpy(line[numberLines++],"0#  pragma warning(disable:4786)");
4563  strcpy(line[numberLines++],"0#endif\n");
4564  strcpy(line[numberLines++],"0#include <cassert>");
4565  strcpy(line[numberLines++],"0#include <iomanip>");
4566  strcpy(line[numberLines++],"0#include \"OsiClpSolverInterface.hpp\"");
4567  strcpy(line[numberLines++],"0#include \"CbcModel.hpp\"");
4568  strcpy(line[numberLines++],"0#include \"CbcCutGenerator.hpp\"");
4569  strcpy(line[numberLines++],"0#include \"CbcStrategy.hpp\"");
4570  strcpy(line[numberLines++],"0#include \"CglPreProcess.hpp\"");
4571  strcpy(line[numberLines++],"0#include \"CoinTime.hpp\"");
4572  if (preProcess>0) 
4573    strcpy(line[numberLines++],"0#include \"CglProbing.hpp\""); // possibly redundant
4574  // To allow generated 5's to be just before branchAndBound - do rest here
4575  strcpy(line[numberLines++],"5  cbcModel->initialSolve();");
4576  strcpy(line[numberLines++],"5  if (clpModel->tightenPrimalBounds()!=0) {");
4577  strcpy(line[numberLines++],"5    std::cout<<\"Problem is infeasible - tightenPrimalBounds!\"<<std::endl;");
4578  strcpy(line[numberLines++],"5    exit(1);");
4579  strcpy(line[numberLines++],"5  }");
4580  strcpy(line[numberLines++],"5  clpModel->dual();  // clean up");
4581  if (sizecode) {
4582    // override some settings
4583    strcpy(line[numberLines++],"5  // compute some things using problem size");
4584    strcpy(line[numberLines++],"5  cbcModel->setMinimumDrop(min(5.0e-2,");
4585    strcpy(line[numberLines++],"5       fabs(cbcModel->getMinimizationObjValue())*1.0e-3+1.0e-4));");
4586    strcpy(line[numberLines++],"5  if (cbcModel->getNumCols()<500)");
4587    strcpy(line[numberLines++],"5    cbcModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible");
4588    strcpy(line[numberLines++],"5  else if (cbcModel->getNumCols()<5000)");
4589    strcpy(line[numberLines++],"5    cbcModel->setMaximumCutPassesAtRoot(100); // use minimum drop");
4590    strcpy(line[numberLines++],"5  else");
4591    strcpy(line[numberLines++],"5    cbcModel->setMaximumCutPassesAtRoot(20);");
4592    strcpy(line[numberLines++],"5  cbcModel->setMaximumCutPasses(1);");
4593  }
4594  if (preProcess<=0) {
4595    // no preprocessing or strategy
4596    if (preProcess) {
4597      strcpy(line[numberLines++],"5  // Preprocessing using CbcStrategy");
4598      strcpy(line[numberLines++],"5  CbcStrategyDefault strategy(true,5,5);");
4599      strcpy(line[numberLines++],"5  strategy.setupPreProcessing(1);");
4600      strcpy(line[numberLines++],"5  cbcModel->setStrategy(strategy);");
4601    }
4602  } else {
4603    int translate[]={9999,0,0,-1,2,3,-2};
4604    strcpy(line[numberLines++],"5  // Hand coded preprocessing");
4605    strcpy(line[numberLines++],"5  CglPreProcess process;");
4606    strcpy(line[numberLines++],"5  OsiSolverInterface * saveSolver=cbcModel->solver()->clone();");
4607    strcpy(line[numberLines++],"5  // Tell solver we are in Branch and Cut");
4608    strcpy(line[numberLines++],"5  saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;");
4609    strcpy(line[numberLines++],"5  // Default set of cut generators");
4610    strcpy(line[numberLines++],"5  CglProbing generator1;");
4611    strcpy(line[numberLines++],"5  generator1.setUsingObjective(true);");
4612    strcpy(line[numberLines++],"5  generator1.setMaxPass(3);");
4613    strcpy(line[numberLines++],"5  generator1.setMaxProbeRoot(saveSolver->getNumCols());");
4614    strcpy(line[numberLines++],"5  generator1.setMaxElements(100);");
4615    strcpy(line[numberLines++],"5  generator1.setMaxLookRoot(50);");
4616    strcpy(line[numberLines++],"5  generator1.setRowCuts(3);");
4617    strcpy(line[numberLines++],"5  // Add in generators");
4618    strcpy(line[numberLines++],"5  process.addCutGenerator(&generator1);");
4619    strcpy(line[numberLines++],"5  process.messageHandler()->setLogLevel(cbcModel->logLevel());");
4620    strcpy(line[numberLines++],"5  OsiSolverInterface * solver2 = ");
4621    sprintf(line[numberLines++],"5    process.preProcessNonDefault(*saveSolver,%d,10);",translate[preProcess]);
4622    strcpy(line[numberLines++],"5  // Tell solver we are not in Branch and Cut");
4623    strcpy(line[numberLines++],"5  saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;");
4624    strcpy(line[numberLines++],"5  if (solver2)");
4625    strcpy(line[numberLines++],"5    solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;");
4626    strcpy(line[numberLines++],"5  if (!solver2) {");
4627    strcpy(line[numberLines++],"5    std::cout<<\"Pre-processing says infeasible!\"<<std::endl;");
4628    strcpy(line[numberLines++],"5    exit(1);");
4629    strcpy(line[numberLines++],"5  } else {");
4630    strcpy(line[numberLines++],"5    std::cout<<\"processed model has \"<<solver2->getNumRows()");
4631    strcpy(line[numberLines++],"5            <<\" rows, \"<<solver2->getNumCols()");
4632    strcpy(line[numberLines++],"5            <<\" columns and \"<<solver2->getNumElements()");
4633    strcpy(line[numberLines++],"5            <<\" elements\"<<solver2->getNumElements()<<std::endl;");
4634    strcpy(line[numberLines++],"5  }");
4635    strcpy(line[numberLines++],"5  // we have to keep solver2 so pass clone");
4636    strcpy(line[numberLines++],"5  solver2 = solver2->clone();");
4637    strcpy(line[numberLines++],"5  cbcModel->assignSolver(solver2);");
4638    strcpy(line[numberLines++],"5  cbcModel->initialSolve();");
4639  }
4640  while (fgets(line[numberLines],MAXONELINE,fp)) {
4641    assert (numberLines<MAXLINES);
4642    clean(line[numberLines]);
4643    numberLines++;
4644  }
4645  fclose(fp);
4646  strcpy(line[numberLines++],"0\nint main (int argc, const char *argv[])\n{");
4647  strcpy(line[numberLines++],"0  OsiClpSolverInterface solver1;");
4648  strcpy(line[numberLines++],"0  int status=1;");
4649  strcpy(line[numberLines++],"0  if (argc<2)");
4650  strcpy(line[numberLines++],"0    std::cout<<\"Please give file name\"<<std::endl;");
4651  strcpy(line[numberLines++],"0  else");
4652  strcpy(line[numberLines++],"0    status=solver1.readMps(argv[1],\"\");");
4653  strcpy(line[numberLines++],"0  if (status) {");
4654  strcpy(line[numberLines++],"0    std::cout<<\"Bad readMps \"<<argv[1]<<std::endl;");
4655  strcpy(line[numberLines++],"0    exit(1);");
4656  strcpy(line[numberLines++],"0  }\n");
4657  strcpy(line[numberLines++],"0  double time1 = CoinCpuTime();");
4658  strcpy(line[numberLines++],"0  CbcModel model(solver1);");
4659  strcpy(line[numberLines++],"0  // Now do requested saves and modifications");
4660  strcpy(line[numberLines++],"0  CbcModel * cbcModel = & model;");
4661  strcpy(line[numberLines++],"0  OsiSolverInterface * osiModel = model.solver();");
4662  strcpy(line[numberLines++],"0  OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);");
4663  strcpy(line[numberLines++],"0  ClpSimplex * clpModel = osiclpModel->getModelPtr();");
4664  // add in comments about messages
4665  strcpy(line[numberLines++],"3  // You can save some time by switching off message building");
4666  strcpy(line[numberLines++],"3  // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);");
4667  // add in actual solve
4668  strcpy(line[numberLines++],"5  cbcModel->branchAndBound();");
4669  strcpy(line[numberLines++],"8  std::cout<<argv[1]<<\" took \"<<CoinCpuTime()-time1<<\" seconds, \"");
4670  strcpy(line[numberLines++],"8    <<cbcModel->getNodeCount()<<\" nodes with objective \"");
4671  strcpy(line[numberLines++],"8    <<cbcModel->getObjValue()");
4672  strcpy(line[numberLines++],"8    <<(!cbcModel->status() ? \" Finished\" : \" Not finished\")");
4673  strcpy(line[numberLines++],"8    <<std::endl;");
4674  strcpy(line[numberLines++],"5  // For best solution");
4675  strcpy(line[numberLines++],"5  int numberColumns = solver1.getNumCols();");
4676  strcpy(line[numberLines++],"5  if (cbcModel->getMinimizationObjValue()<1.0e50) {");
4677  if (preProcess>0) {
4678    strcpy(line[numberLines++],"5    // post process");
4679    strcpy(line[numberLines++],"5    process.postProcess(*cbcModel->solver());");
4680    strcpy(line[numberLines++],"5    // Solution now back in saveSolver");
4681    strcpy(line[numberLines++],"5    cbcModel->assignSolver(saveSolver);");
4682    strcpy(line[numberLines++],"5    memcpy(cbcModel->bestSolution(),cbcModel->solver()->getColSolution(),");
4683    strcpy(line[numberLines++],"5          numberColumns*sizeof(double));");
4684  }
4685  strcpy(line[numberLines++],"5    // put back in original solver");
4686  strcpy(line[numberLines++],"5    solver1.setColSolution(cbcModel->bestSolution());");
4687  strcpy(line[numberLines++],"5    const double * solution = solver1.getColSolution();");
4688  strcpy(line[numberLines++],"8  \n  // Now you would use solution etc etc\n");
4689  strcpy(line[numberLines++],"5");
4690  strcpy(line[numberLines++],"5    // Get names from solver1 (as OsiSolverInterface may lose)");
4691  strcpy(line[numberLines++],"5    std::vector<std::string> columnNames = *solver1.getModelPtr()->columnNames();");
4692  strcpy(line[numberLines++],"5    ");
4693  strcpy(line[numberLines++],"5    int iColumn;");
4694  strcpy(line[numberLines++],"5    std::cout<<std::setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14);");
4695  strcpy(line[numberLines++],"5    ");
4696  strcpy(line[numberLines++],"5    std::cout<<\"--------------------------------------\"<<std::endl;");
4697  strcpy(line[numberLines++],"5    for (iColumn=0;iColumn<numberColumns;iColumn++) {");
4698  strcpy(line[numberLines++],"5      double value=solution[iColumn];");
4699  strcpy(line[numberLines++],"5      if (fabs(value)>1.0e-7&&solver1.isInteger(iColumn)) ");
4700  strcpy(line[numberLines++],"5 std::cout<<std::setw(6)<<iColumn<<\" \"");
4701  strcpy(line[numberLines++],"5                 <<columnNames[iColumn]<<\" \"");
4702  strcpy(line[numberLines++],"5                 <<value<<std::endl;");
4703  strcpy(line[numberLines++],"5    }");
4704  strcpy(line[numberLines++],"5    std::cout<<\"--------------------------------------\"<<std::endl;");
4705  strcpy(line[numberLines++],"5  ");
4706  strcpy(line[numberLines++],"5    std::cout<<std::resetiosflags(std::ios::fixed|std::ios::showpoint|std::ios::scientific);");
4707  strcpy(line[numberLines++],"5  }");
4708  strcpy(line[numberLines++],"8  return 0;\n}");
4709  fp = fopen(fileName,"w");
4710  assert (fp);
4711
4712  int wanted[9];
4713  memset(wanted,0,sizeof(wanted));
4714  wanted[0]=wanted[3]=wanted[5]=wanted[8]=1;
4715  if (type>0) 
4716    wanted[1]=wanted[6]=1;
4717  if (type>1) 
4718    wanted[2]=wanted[4]=wanted[7]=1;
4719  std::string header[9]=
4720  { "","Save values","Redundant save of default values","Set changed values",
4721    "Redundant set default values","Solve","Restore values","Redundant restore values","Finish up"};
4722  for (int iType=0;iType<9;iType++) {
4723    if (!wanted[iType])
4724      continue;
4725    int n=0;
4726    int iLine;
4727    for (iLine=0;iLine<numberLines;iLine++) {
4728      if (line[iLine][0]=='0'+iType) {
4729        if (!n&&header[iType]!="")
4730          fprintf(fp,"\n  // %s\n\n",header[iType].c_str());
4731        n++;
4732        // skip save and clp as cloned
4733        if (!strstr(line[iLine],"save")||(!strstr(line[iLine],"clpMo")&&
4734                                          !strstr(line[iLine],"_Osi")))
4735          fprintf(fp,"%s\n",line[iLine]+1);
4736      }
4737    }
4738  }
4739  fclose(fp);
4740  printf("C++ file written to %s\n",fileName);
4741}
4742/*
4743  Version 1.00.00 November 16 2005.
4744  This is to stop me (JJF) messing about too much.
4745  Tuning changes should be noted here.
4746  The testing next version may be activated by CBC_NEXT_VERSION
4747  This applies to OsiClp, Clp etc
4748  Version 1.00.01 November 24 2005
4749  Added several classes for advanced users.  This can't affect code (if you don't use it)
4750  Made some tiny changes (for N way branching) which should not change anything.
4751  CbcNWay object class - for N way branching this also allows use of CbcConsequence class.
4752  CbcBranchAllDifferent object class - for branching on general integer variables
4753  to stop them having same value so branches are x >= y+1 and x <= y-1.
4754  Added two new Cgl classes - CglAllDifferent which does column fixing (too slowly)
4755  and CglStored which just has a list of cuts which can be activated.
4756  Modified preprocess option to SOS
4757  Version 1.00.02 December 9 2005
4758  Added use of CbcStrategy to do clean preprocessing
4759  Added use of referenceSolver for cleaner repetition of Cbc
4760  Version 1.01.00 February 2 2006
4761  Added first try at Ampl interface
4762*/
Note: See TracBrowser for help on using the repository browser.