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

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

for osi methods

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