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

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

extend fpump heuristic and allow more probing

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