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

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

many changes

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