source: trunk/Test/CoinSolve.cpp @ 202

Last change on this file since 202 was 202, checked in by forrest, 16 years ago

try and improve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 83.5 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3   
4#include "CoinPragma.hpp"
5
6#include <cassert>
7#include <cstdio>
8#include <cmath>
9#include <cfloat>
10#include <string>
11#include <iostream>
12
13
14#include "CoinPragma.hpp"
15#include "CoinHelperFunctions.hpp"
16// Same version as CBC
17#define CBCVERSION "0.99"
18
19#include "CoinMpsIO.hpp"
20
21#include "ClpFactorization.hpp"
22#include "CoinTime.hpp"
23#include "ClpSimplex.hpp"
24#include "ClpSolve.hpp"
25#include "ClpPackedMatrix.hpp"
26#include "ClpPlusMinusOneMatrix.hpp"
27#include "ClpNetworkMatrix.hpp"
28#include "ClpDualRowSteepest.hpp"
29#include "ClpDualRowDantzig.hpp"
30#include "ClpLinearObjective.hpp"
31#include "ClpPrimalColumnSteepest.hpp"
32#include "ClpPrimalColumnDantzig.hpp"
33#include "ClpPresolve.hpp"
34#include "CbcOrClpParam.hpp"
35#include "OsiRowCutDebugger.hpp"
36#ifdef DMALLOC
37#include "dmalloc.h"
38#endif
39#ifdef WSSMP_BARRIER
40#define FOREIGN_BARRIER
41#endif
42#ifdef UFL_BARRIER
43#define FOREIGN_BARRIER
44#endif
45#ifdef TAUCS_BARRIER
46#define FOREIGN_BARRIER
47#endif
48#include "CoinWarmStartBasis.hpp"
49
50#include "OsiSolverInterface.hpp"
51#include "OsiCuts.hpp"
52#include "OsiRowCut.hpp"
53#include "OsiColCut.hpp"
54
55#include "CglPreProcess.hpp"
56#include "CglCutGenerator.hpp"
57#include "CglGomory.hpp"
58#include "CglProbing.hpp"
59#include "CglKnapsackCover.hpp"
60#include "CglRedSplit.hpp"
61#include "CglClique.hpp"
62#include "CglFlowCover.hpp"
63#include "CglMixedIntegerRounding2.hpp"
64#include "CglTwomir.hpp"
65
66#include "CbcModel.hpp"
67#include "CbcHeuristic.hpp"
68#include "CbcHeuristicLocal.hpp"
69#include "CbcHeuristicGreedy.hpp"
70#include "CbcHeuristicFPump.hpp"
71#include "CbcTreeLocal.hpp"
72#include "CbcCompareActual.hpp"
73#include  "CbcOrClpParam.hpp"
74#include  "CbcCutGenerator.hpp"
75
76#include "OsiClpSolverInterface.hpp"
77
78static double totalTime=0.0;
79
80//#############################################################################
81
82#ifdef NDEBUG
83#undef NDEBUG
84#endif
85// Allow for interrupts
86// But is this threadsafe ? (so switched off by option)
87
88#include "CoinSignal.hpp"
89static CbcModel * currentBranchModel = NULL;
90
91extern "C" {
92   static void signal_handler(int whichSignal)
93   {
94      if (currentBranchModel!=NULL) 
95         currentBranchModel->setMaximumNodes(0); // stop at next node
96      return;
97   }
98}
99
100int mainTest (int argc, const char *argv[],int algorithm,
101              ClpSimplex empty, bool doPresolve,int doIdiot);
102int CbcOrClpRead_mode=1;
103FILE * CbcOrClpReadCommand=stdin;
104static int * analyze(OsiClpSolverInterface * solverMod, int & numberChanged, double & increment,
105                     bool changeInt)
106{
107  OsiSolverInterface * solver = solverMod->clone();
108  {
109    // just get increment
110    CbcModel model(*solver);
111    model.analyzeObjective();
112    double increment2=model.getCutoffIncrement();
113    printf("initial cutoff increment %g\n",increment2);
114  }
115  const double *objective = solver->getObjCoefficients() ;
116  const double *lower = solver->getColLower() ;
117  const double *upper = solver->getColUpper() ;
118  int numberColumns = solver->getNumCols() ;
119  int numberRows = solver->getNumRows();
120  double direction = solver->getObjSense();
121  int iRow,iColumn;
122
123  // Row copy
124  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
125  const double * elementByRow = matrixByRow.getElements();
126  const int * column = matrixByRow.getIndices();
127  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
128  const int * rowLength = matrixByRow.getVectorLengths();
129
130  // Column copy
131  CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
132  const double * element = matrixByCol.getElements();
133  const int * row = matrixByCol.getIndices();
134  const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
135  const int * columnLength = matrixByCol.getVectorLengths();
136
137  const double * rowLower = solver->getRowLower();
138  const double * rowUpper = solver->getRowUpper();
139
140  char * ignore = new char [numberRows];
141  int * changed = new int[numberColumns];
142  int * which = new int[numberRows];
143  double * changeRhs = new double[numberRows];
144  memset(changeRhs,0,numberRows*sizeof(double));
145  memset(ignore,0,numberRows);
146  numberChanged=0;
147  int numberInteger=0;
148  for (iColumn=0;iColumn<numberColumns;iColumn++) {
149    if (upper[iColumn] > lower[iColumn]+1.0e-8&&solver->isInteger(iColumn)) 
150      numberInteger++;
151  }
152  bool finished=false;
153  while (!finished) {
154    int saveNumberChanged = numberChanged;
155    for (iRow=0;iRow<numberRows;iRow++) {
156      int numberContinuous=0;
157      double value1=0.0,value2=0.0;
158      bool allIntegerCoeff=true;
159      double sumFixed=0.0;
160      int jColumn1=-1,jColumn2=-1;
161      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
162        int jColumn = column[j];
163        double value = elementByRow[j];
164        if (upper[jColumn] > lower[jColumn]+1.0e-8) {
165          if (!solver->isInteger(jColumn)) {
166            if (numberContinuous==0) {
167              jColumn1=jColumn;
168              value1=value;
169            } else {
170              jColumn2=jColumn;
171              value2=value;
172            }
173            numberContinuous++;
174          } else {
175            if (fabs(value-floor(value+0.5))>1.0e-12)
176              allIntegerCoeff=false;
177          }
178        } else {
179          sumFixed += lower[jColumn]*value;
180        }
181      }
182      double low = rowLower[iRow];
183      if (low>-1.0e20) {
184        low -= sumFixed;
185        if (fabs(low-floor(low+0.5))>1.0e-12)
186          allIntegerCoeff=false;
187      }
188      double up = rowUpper[iRow];
189      if (up<1.0e20) {
190        up -= sumFixed;
191        if (fabs(up-floor(up+0.5))>1.0e-12)
192          allIntegerCoeff=false;
193      }
194      if (numberContinuous==1) {
195        // see if really integer
196        // This does not allow for complicated cases
197        if (low==up) {
198          if (fabs(value1)>1.0e-3) {
199            value1 = 1.0/value1;
200            if (fabs(value1-floor(value1+0.5))<1.0e-12) {
201              // integer
202              changed[numberChanged++]=jColumn1;
203              solver->setInteger(jColumn1);
204              if (upper[jColumn1]>1.0e20)
205                solver->setColUpper(jColumn1,1.0e20);
206              if (lower[jColumn1]<-1.0e20)
207                solver->setColLower(jColumn1,-1.0e20);
208            }
209          }
210        } else {
211          if (fabs(value1)>1.0e-3) {
212            value1 = 1.0/value1;
213            if (fabs(value1-floor(value1+0.5))<1.0e-12) {
214              // This constraint will not stop it being integer
215              ignore[iRow]=1;
216            }
217          }
218        }
219      } else if (numberContinuous==2) {
220        if (low==up) {
221          /* need general theory - for now just look at 2 cases -
222             1 - +- 1 one in column and just costs i.e. matching objective
223             2 - +- 1 two in column but feeds into G/L row which will try and minimize
224          */
225          if (fabs(value1)==1.0&&value1*value2==-1.0&&!lower[jColumn1]
226              &&!lower[jColumn2]) {
227            int n=0;
228            int i;
229            double objChange=direction*(objective[jColumn1]+objective[jColumn2]);
230            double bound = CoinMin(upper[jColumn1],upper[jColumn2]);
231            bound = CoinMin(bound,1.0e20);
232            for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
233              int jRow = row[i];
234              double value = element[i];
235              if (jRow!=iRow) {
236                which[n++]=jRow;
237                changeRhs[jRow]=value;
238              }
239            }
240            for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
241              int jRow = row[i];
242              double value = element[i];
243              if (jRow!=iRow) {
244                if (!changeRhs[jRow]) {
245                  which[n++]=jRow;
246                  changeRhs[jRow]=value;
247                } else {
248                  changeRhs[jRow]+=value;
249                }
250              }
251            }
252            if (objChange>=0.0) {
253              // see if all rows OK
254              bool good=true;
255              for (i=0;i<n;i++) {
256                int jRow = which[i];
257                double value = changeRhs[jRow];
258                if (value) {
259                  value *= bound;
260                  if (rowLength[jRow]==1) {
261                    if (value>0.0) {
262                      double rhs = rowLower[jRow];
263                      if (rhs>0.0) {
264                        double ratio =rhs/value;
265                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
266                          good=false;
267                      }
268                    } else {
269                      double rhs = rowUpper[jRow];
270                      if (rhs<0.0) {
271                        double ratio =rhs/value;
272                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
273                          good=false;
274                      }
275                    }
276                  } else if (rowLength[jRow]==2) {
277                    if (value>0.0) {
278                      if (rowLower[jRow]>-1.0e20)
279                        good=false;
280                    } else {
281                      if (rowUpper[jRow]<1.0e20)
282                        good=false;
283                    }
284                  } else {
285                    good=false;
286                  }
287                }
288              }
289              if (good) {
290                // both can be integer
291                changed[numberChanged++]=jColumn1;
292                solver->setInteger(jColumn1);
293                if (upper[jColumn1]>1.0e20)
294                  solver->setColUpper(jColumn1,1.0e20);
295                if (lower[jColumn1]<-1.0e20)
296                  solver->setColLower(jColumn1,-1.0e20);
297                changed[numberChanged++]=jColumn2;
298                solver->setInteger(jColumn2);
299                if (upper[jColumn2]>1.0e20)
300                  solver->setColUpper(jColumn2,1.0e20);
301                if (lower[jColumn2]<-1.0e20)
302                  solver->setColLower(jColumn2,-1.0e20);
303              }
304            }
305            // clear
306            for (i=0;i<n;i++) {
307              changeRhs[which[i]]=0.0;
308            }
309          }
310        }
311      }
312    }
313    for (iColumn=0;iColumn<numberColumns;iColumn++) {
314      if (upper[iColumn] > lower[iColumn]+1.0e-8&&!solver->isInteger(iColumn)) {
315        double value;
316        value = upper[iColumn];
317        if (value<1.0e20&&fabs(value-floor(value+0.5))>1.0e-12) 
318          continue;
319        value = lower[iColumn];
320        if (value>-1.0e20&&fabs(value-floor(value+0.5))>1.0e-12) 
321          continue;
322        bool integer=true;
323        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
324          int iRow = row[j];
325          if (!ignore[iRow]) {
326            integer=false;
327            break;
328          }
329        }
330        if (integer) {
331          // integer
332          changed[numberChanged++]=iColumn;
333          solver->setInteger(iColumn);
334          if (upper[iColumn]>1.0e20)
335            solver->setColUpper(iColumn,1.0e20);
336          if (lower[iColumn]<-1.0e20)
337            solver->setColLower(iColumn,-1.0e20);
338        }
339      }
340    }
341    finished = numberChanged==saveNumberChanged;
342  }
343  delete [] which;
344  delete [] changeRhs;
345  if (numberInteger)
346    printf("%d integer variables",numberInteger);
347  if (changeInt) {
348    if (numberChanged)
349      printf(" and %d variables made integer\n",numberChanged);
350    else
351      printf("\n");
352    delete [] ignore;
353    //increment=0.0;
354    if (!numberChanged) {
355      delete [] changed;
356      delete solver;
357      return NULL;
358    } else {
359      for (iColumn=0;iColumn<numberColumns;iColumn++) {
360        if (solver->isInteger(iColumn))
361          solverMod->setInteger(iColumn);
362      }
363      delete solver;
364      return changed;
365    }
366  } else {
367    if (numberChanged)
368      printf(" and %d variables could be made integer\n",numberChanged);
369    else
370      printf("\n");
371    // just get increment
372    CbcModel model(*solver);
373    model.analyzeObjective();
374    double increment2=model.getCutoffIncrement();
375    if (increment2>increment) {
376      printf("cutoff increment increased from %g to %g\n",increment,increment2);
377      increment=increment2;
378    }
379    delete solver;
380    numberChanged=0;
381    return NULL;
382  }
383}
384int main (int argc, const char *argv[])
385{
386  /* Note
387     This is meant as a stand-alone executable to do as much of coin as possible.
388     It should only have one solver known to it.
389  */
390  {
391    double time1 = CoinCpuTime(),time2;
392    CoinSighandler_t saveSignal=SIG_DFL;
393    // register signal handler
394    saveSignal = signal(SIGINT,signal_handler);
395    // Set up all non-standard stuff
396    OsiClpSolverInterface solver1;
397    CbcModel model(solver1);
398    CbcModel * babModel = NULL;
399    model.setNumberBeforeTrust(21);
400    int cutPass=-1234567;
401    OsiSolverInterface * solver = model.solver();
402    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
403    ClpSimplex * lpSolver = clpSolver->getModelPtr();
404    clpSolver->messageHandler()->setLogLevel(0) ;
405    model.messageHandler()->setLogLevel(1);
406   
407   
408    // default action on import
409    int allowImportErrors=0;
410    int keepImportNames=1;
411    int doIdiot=-1;
412    int outputFormat=2;
413    int slpValue=-1;
414    int printOptions=0;
415    int printMode=0;
416    int presolveOptions=0;
417    int doCrash=0;
418    int doSprint=-1;
419    int doScaling=1;
420    // set reasonable defaults
421    int preSolve=5;
422    int preProcess=1;
423    bool preSolveFile=false;
424   
425    double djFix=1.0e100;
426    double gapRatio=1.0e100;
427    double tightenFactor=0.0;
428    lpSolver->setPerturbation(50);
429    lpSolver->messageHandler()->setPrefix(false);
430    const char dirsep =  CoinFindDirSeparator();
431    std::string directory = (dirsep == '/' ? "./" : ".\\");
432    std::string defaultDirectory = directory;
433    std::string importFile ="";
434    std::string exportFile ="default.mps";
435    std::string importBasisFile ="";
436    std::string debugFile="";
437    double * debugValues = NULL;
438    int numberDebugValues = -1;
439    int basisHasValues=0;
440    std::string exportBasisFile ="default.bas";
441    std::string saveFile ="default.prob";
442    std::string restoreFile ="default.prob";
443    std::string solutionFile ="stdout";
444    std::string solutionSaveFile ="solution.file";
445#define CBCMAXPARAMETERS 200
446    CbcOrClpParam parameters[CBCMAXPARAMETERS];
447    int numberParameters ;
448    establishParams(numberParameters,parameters) ;
449    parameters[whichParam(BASISIN,numberParameters,parameters)].setStringValue(importBasisFile);
450    parameters[whichParam(BASISOUT,numberParameters,parameters)].setStringValue(exportBasisFile);
451    parameters[whichParam(DEBUG,numberParameters,parameters)].setStringValue(debugFile);
452    parameters[whichParam(DIRECTORY,numberParameters,parameters)].setStringValue(directory);
453    parameters[whichParam(DUALBOUND,numberParameters,parameters)].setDoubleValue(lpSolver->dualBound());
454    parameters[whichParam(DUALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->dualTolerance());
455    parameters[whichParam(EXPORT,numberParameters,parameters)].setStringValue(exportFile);
456    parameters[whichParam(IDIOT,numberParameters,parameters)].setIntValue(doIdiot);
457    parameters[whichParam(IMPORT,numberParameters,parameters)].setStringValue(importFile);
458    int slog = whichParam(SOLVERLOGLEVEL,numberParameters,parameters);
459    int log = whichParam(LOGLEVEL,numberParameters,parameters);
460    parameters[slog].setIntValue(0);
461    parameters[log].setIntValue(1);
462    parameters[whichParam(MAXFACTOR,numberParameters,parameters)].setIntValue(lpSolver->factorizationFrequency());
463    parameters[whichParam(MAXITERATION,numberParameters,parameters)].setIntValue(lpSolver->maximumIterations());
464    parameters[whichParam(OUTPUTFORMAT,numberParameters,parameters)].setIntValue(outputFormat);
465    parameters[whichParam(PRESOLVEPASS,numberParameters,parameters)].setIntValue(preSolve);
466    parameters[whichParam(PERTVALUE,numberParameters,parameters)].setIntValue(lpSolver->perturbation());
467    parameters[whichParam(PRIMALTOLERANCE,numberParameters,parameters)].setDoubleValue(lpSolver->primalTolerance());
468    parameters[whichParam(PRIMALWEIGHT,numberParameters,parameters)].setDoubleValue(lpSolver->infeasibilityCost());
469    parameters[whichParam(RESTORE,numberParameters,parameters)].setStringValue(restoreFile);
470    parameters[whichParam(SAVE,numberParameters,parameters)].setStringValue(saveFile);
471    //parameters[whichParam(TIMELIMIT,numberParameters,parameters)].setDoubleValue(1.0e8);
472    parameters[whichParam(TIMELIMIT_BAB,numberParameters,parameters)].setDoubleValue(1.0e8);
473    parameters[whichParam(SOLUTION,numberParameters,parameters)].setStringValue(solutionFile);
474    parameters[whichParam(SAVESOL,numberParameters,parameters)].setStringValue(solutionSaveFile);
475    parameters[whichParam(SPRINT,numberParameters,parameters)].setIntValue(doSprint);
476    model.setNumberBeforeTrust(5);
477    parameters[whichParam(NUMBERBEFORE,numberParameters,parameters)].setIntValue(5);
478    parameters[whichParam(MAXNODES,numberParameters,parameters)].setIntValue(model.getMaximumNodes());
479    model.setNumberStrong(5);
480    parameters[whichParam(STRONGBRANCHING,numberParameters,parameters)].setIntValue(model.numberStrong());
481    parameters[whichParam(INFEASIBILITYWEIGHT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcInfeasibilityWeight));
482    parameters[whichParam(INTEGERTOLERANCE,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcIntegerTolerance));
483    parameters[whichParam(INCREMENT,numberParameters,parameters)].setDoubleValue(model.getDblParam(CbcModel::CbcCutoffIncrement));
484    // Set up likely cut generators and defaults
485    parameters[whichParam(PREPROCESS,numberParameters,parameters)].setCurrentOption("on");
486
487    CglGomory gomoryGen;
488    // try larger limit
489    gomoryGen.setLimitAtRoot(500);
490    gomoryGen.setLimit(100);
491    // set default action (0=off,1=on,2=root)
492    int gomoryAction=3;
493    parameters[whichParam(GOMORYCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
494
495    CglProbing probingGen;
496    probingGen.setUsingObjective(true);
497    probingGen.setMaxPass(3);
498    probingGen.setMaxPassRoot(3);
499    // Number of unsatisfied variables to look at
500    probingGen.setMaxProbe(10);
501    probingGen.setMaxProbeRoot(50);
502    // How far to follow the consequences
503    probingGen.setMaxLook(10);
504    probingGen.setMaxLookRoot(50);
505    probingGen.setMaxLookRoot(10);
506    // Only look at rows with fewer than this number of elements
507    probingGen.setMaxElements(200);
508    probingGen.setRowCuts(3);
509    // set default action (0=off,1=on,2=root)
510    int probingAction=3;
511    parameters[whichParam(PROBINGCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
512
513    CglKnapsackCover knapsackGen;
514    // set default action (0=off,1=on,2=root)
515    int knapsackAction=3;
516    parameters[whichParam(KNAPSACKCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
517
518    CglRedSplit redsplitGen;
519    redsplitGen.setLimit(100);
520    // set default action (0=off,1=on,2=root)
521    // Off as seems to give some bad cuts
522    int redsplitAction=2;
523    parameters[whichParam(REDSPLITCUTS,numberParameters,parameters)].setCurrentOption("root");
524
525    CglClique cliqueGen(false,true);
526    cliqueGen.setStarCliqueReport(false);
527    cliqueGen.setRowCliqueReport(false);
528    cliqueGen.setMinViolation(0.1);
529    // set default action (0=off,1=on,2=root)
530    int cliqueAction=3;
531    parameters[whichParam(CLIQUECUTS,numberParameters,parameters)].setCurrentOption("ifmove");
532
533    CglMixedIntegerRounding2 mixedGen;
534    // set default action (0=off,1=on,2=root)
535    int mixedAction=3;
536    parameters[whichParam(MIXEDCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
537
538    CglFlowCover flowGen;
539    // set default action (0=off,1=on,2=root)
540    int flowAction=3;
541    parameters[whichParam(FLOWCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
542
543    CglTwomir twomirGen;
544    // set default action (0=off,1=on,2=root)
545    int twomirAction=3;
546    parameters[whichParam(TWOMIRCUTS,numberParameters,parameters)].setCurrentOption("ifmove");
547
548    bool useRounding=true;
549    parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption("on");
550    bool useFpump=true;
551    bool useGreedy=true;
552    parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption("on");
553    bool useCombine=true;
554    parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption("on");
555    bool useLocalTree=false;
556   
557    // total number of commands read
558    int numberGoodCommands=0;
559    bool goodModel=false;
560    // Set false if user does anything advanced
561    bool defaultSettings=true;
562
563    // Hidden stuff for barrier
564    int choleskyType = 0;
565    int gamma=0;
566    int scaleBarrier=0;
567    int doKKT=0;
568    int crossover=2; // do crossover unless quadratic
569    // For names
570    int lengthName = 0;
571    std::vector<std::string> rowNames;
572    std::vector<std::string> columnNames;
573   
574    std::string field;
575    std::cout<<"Coin Cbc and Clp Solver version "<<CBCVERSION
576             <<", build "<<__DATE__<<std::endl;
577   
578    while (1) {
579      // next command
580      field=CoinReadGetCommand(argc,argv);
581     
582      // exit if null or similar
583      if (!field.length()) {
584        if (numberGoodCommands==1&&goodModel) {
585          // we just had file name - do branch and bound
586          field="branch";
587        } else if (!numberGoodCommands) {
588          // let's give the sucker a hint
589          std::cout
590            <<"CoinSolver takes input from arguments ( - switches to stdin)"
591            <<std::endl
592            <<"Enter ? for list of commands or help"<<std::endl;
593          field="-";
594        } else {
595          break;
596        }
597      }
598     
599      // see if ? at end
600      int numberQuery=0;
601      if (field!="?"&&field!="???") {
602        int length = field.length();
603        int i;
604        for (i=length-1;i>0;i--) {
605          if (field[i]=='?') 
606            numberQuery++;
607          else
608            break;
609        }
610        field=field.substr(0,length-numberQuery);
611      }
612      // find out if valid command
613      int iParam;
614      int numberMatches=0;
615      int firstMatch=-1;
616      for ( iParam=0; iParam<numberParameters; iParam++ ) {
617        int match = parameters[iParam].matches(field);
618        if (match==1) {
619          numberMatches = 1;
620          firstMatch=iParam;
621          break;
622        } else {
623          if (match&&firstMatch<0)
624            firstMatch=iParam;
625          numberMatches += match>>1;
626        }
627      }
628      if (iParam<numberParameters&&!numberQuery) {
629        // found
630        CbcOrClpParam found = parameters[iParam];
631        CbcOrClpParameterType type = found.type();
632        int valid;
633        numberGoodCommands++;
634        if (type==BAB&&goodModel) {
635          // check if any integers
636          if (!lpSolver->integerInformation())
637            type=DUALSIMPLEX;
638        }
639        if (type==GENERALQUERY) {
640          std::cout<<"In argument list keywords have leading - "
641            ", -stdin or just - switches to stdin"<<std::endl;
642          std::cout<<"One command per line (and no -)"<<std::endl;
643          std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
644          std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
645          std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
646          std::cout<<"abcd value sets value"<<std::endl;
647          std::cout<<"Commands are:"<<std::endl;
648          int maxAcross=5;
649          int limits[]={1,51,101,151,201,251,301,351,401};
650          std::vector<std::string> types;
651          types.push_back("Double parameters:");
652          types.push_back("Branch and Cut double parameters:");
653          types.push_back("Integer parameters:");
654          types.push_back("Branch and Cut integer parameters:");
655          types.push_back("Keyword parameters and others:");
656          types.push_back("Branch and Cut keyword parameters and others:");
657          types.push_back("Actions:");
658          types.push_back("Branch and Cut actions:");
659          int iType;
660          for (iType=0;iType<8;iType++) {
661            int across=0;
662            std::cout<<types[iType]<<"   ";
663            for ( iParam=0; iParam<numberParameters; iParam++ ) {
664              int type = parameters[iParam].type();
665              if (parameters[iParam].displayThis()&&type>=limits[iType]
666                  &&type<limits[iType+1]) {
667                if (!across)
668                  std::cout<<"  ";
669                std::cout<<parameters[iParam].matchName()<<"  ";
670                across++;
671                if (across==maxAcross) {
672                  std::cout<<std::endl;
673                  across=0;
674                }
675              }
676            }
677            if (across)
678              std::cout<<std::endl;
679          }
680        } else if (type==FULLGENERALQUERY) {
681          std::cout<<"Full list of commands is:"<<std::endl;
682          int maxAcross=5;
683          int limits[]={1,51,101,151,201,251,301,351,401};
684          std::vector<std::string> types;
685          types.push_back("Double parameters:");
686          types.push_back("Branch and Cut double parameters:");
687          types.push_back("Integer parameters:");
688          types.push_back("Branch and Cut integer parameters:");
689          types.push_back("Keyword parameters and others:");
690          types.push_back("Branch and Cut keyword parameters and others:");
691          types.push_back("Actions:");
692          types.push_back("Branch and Cut actions:");
693          int iType;
694          for (iType=0;iType<8;iType++) {
695            int across=0;
696            std::cout<<types[iType]<<"  ";
697            for ( iParam=0; iParam<numberParameters; iParam++ ) {
698              int type = parameters[iParam].type();
699              if (type>=limits[iType]
700                  &&type<limits[iType+1]) {
701                if (!across)
702                  std::cout<<"  ";
703                std::cout<<parameters[iParam].matchName()<<"  ";
704                across++;
705                if (across==maxAcross) {
706                  std::cout<<std::endl;
707                  across=0;
708                }
709              }
710            }
711            if (across)
712              std::cout<<std::endl;
713          }
714        } else if (type<101) {
715          // get next field as double
716          double value = CoinReadGetDoubleField(argc,argv,&valid);
717          if (!valid) {
718            parameters[iParam].setDoubleValue(value);
719            if (type<51) {
720              parameters[iParam].setDoubleParameter(lpSolver,value);
721            } else if (type<81) {
722              parameters[iParam].setDoubleParameter(model,value);
723            } else {
724              switch(type) {
725              case DJFIX:
726                djFix=value;
727                preSolve=5;
728                defaultSettings=false; // user knows what she is doing
729                break;
730              case GAPRATIO:
731                gapRatio=value;
732                break;
733              case TIGHTENFACTOR:
734                tightenFactor=value;
735                defaultSettings=false; // user knows what she is doing
736                break;
737              default:
738                abort();
739              }
740            }
741          } else if (valid==1) {
742            abort();
743          } else {
744            std::cout<<parameters[iParam].name()<<" has value "<<
745              parameters[iParam].doubleValue()<<std::endl;
746          }
747        } else if (type<201) {
748          // get next field as int
749          int value = CoinReadGetIntField(argc,argv,&valid);
750          if (!valid) {
751            parameters[iParam].setIntValue(value);
752            if (type<151) {
753              if (parameters[iParam].type()==PRESOLVEPASS)
754                preSolve = value;
755              else if (parameters[iParam].type()==IDIOT)
756                doIdiot = value;
757              else if (parameters[iParam].type()==SPRINT)
758                doSprint = value;
759              else if (parameters[iParam].type()==OUTPUTFORMAT)
760                outputFormat = value;
761              else if (parameters[iParam].type()==SLPVALUE)
762                slpValue = value;
763              else if (parameters[iParam].type()==PRESOLVEOPTIONS)
764                presolveOptions = value;
765              else if (parameters[iParam].type()==PRINTOPTIONS)
766                printOptions = value;
767              else if (parameters[iParam].type()==CUTPASS)
768                cutPass = value;
769              else if (parameters[iParam].type()==FPUMPITS)
770                { useFpump = true;parameters[iParam].setIntValue(value);}
771              else
772                parameters[iParam].setIntParameter(lpSolver,value);
773            } else {
774              parameters[iParam].setIntParameter(model,value);
775            }
776          } else if (valid==1) {
777            abort();
778          } else {
779            std::cout<<parameters[iParam].name()<<" has value "<<
780              parameters[iParam].intValue()<<std::endl;
781          }
782        } else if (type<301) {
783          // one of several strings
784          std::string value = CoinReadGetString(argc,argv);
785          int action = parameters[iParam].parameterOption(value);
786          if (action<0) {
787            if (value!="EOL") {
788              // no match
789              parameters[iParam].printOptions();
790            } else {
791              // print current value
792              std::cout<<parameters[iParam].name()<<" has value "<<
793                parameters[iParam].currentOption()<<std::endl;
794            }
795          } else {
796            parameters[iParam].setCurrentOption(action);
797            // for now hard wired
798            switch (type) {
799            case DIRECTION:
800              if (action==0)
801                lpSolver->setOptimizationDirection(1);
802              else if (action==1)
803                lpSolver->setOptimizationDirection(-1);
804              else
805                lpSolver->setOptimizationDirection(0);
806              break;
807            case DUALPIVOT:
808              if (action==0) {
809                ClpDualRowSteepest steep(3);
810                lpSolver->setDualRowPivotAlgorithm(steep);
811              } else if (action==1) {
812                //ClpDualRowDantzig dantzig;
813                ClpDualRowSteepest dantzig(5);
814                lpSolver->setDualRowPivotAlgorithm(dantzig);
815              } else if (action==2) {
816                // partial steep
817                ClpDualRowSteepest steep(2);
818                lpSolver->setDualRowPivotAlgorithm(steep);
819              } else {
820                ClpDualRowSteepest steep;
821                lpSolver->setDualRowPivotAlgorithm(steep);
822              }
823              break;
824            case PRIMALPIVOT:
825              if (action==0) {
826                ClpPrimalColumnSteepest steep(3);
827                lpSolver->setPrimalColumnPivotAlgorithm(steep);
828              } else if (action==1) {
829                ClpPrimalColumnSteepest steep(0);
830                lpSolver->setPrimalColumnPivotAlgorithm(steep);
831              } else if (action==2) {
832                ClpPrimalColumnDantzig dantzig;
833                lpSolver->setPrimalColumnPivotAlgorithm(dantzig);
834              } else if (action==3) {
835                ClpPrimalColumnSteepest steep(2);
836                lpSolver->setPrimalColumnPivotAlgorithm(steep);
837              } else if (action==4) {
838                ClpPrimalColumnSteepest steep(1);
839                lpSolver->setPrimalColumnPivotAlgorithm(steep);
840              } else if (action==5) {
841                ClpPrimalColumnSteepest steep(4);
842                lpSolver->setPrimalColumnPivotAlgorithm(steep);
843              } else if (action==6) {
844                ClpPrimalColumnSteepest steep(10);
845                lpSolver->setPrimalColumnPivotAlgorithm(steep);
846              }
847              break;
848            case SCALING:
849              lpSolver->scaling(action);
850              solver->setHintParam(OsiDoScale,action!=0,OsiHintTry);
851              doScaling = 1-action;
852              break;
853            case AUTOSCALE:
854              lpSolver->setAutomaticScaling(action!=0);
855              break;
856            case SPARSEFACTOR:
857              lpSolver->setSparseFactorization((1-action)!=0);
858              break;
859            case BIASLU:
860              lpSolver->factorization()->setBiasLU(action);
861              break;
862            case PERTURBATION:
863              if (action==0)
864                lpSolver->setPerturbation(50);
865              else
866                lpSolver->setPerturbation(100);
867              break;
868            case ERRORSALLOWED:
869              allowImportErrors = action;
870              break;
871            case INTPRINT:
872              printMode=action;
873              break;
874              //case ALGORITHM:
875              //algorithm  = action;
876              //defaultSettings=false; // user knows what she is doing
877              //abort();
878              //break;
879            case KEEPNAMES:
880              keepImportNames = 1-action;
881              break;
882            case PRESOLVE:
883              if (action==0)
884                preSolve = 5;
885              else if (action==1)
886                preSolve=0;
887              else if (action==2)
888                preSolve=10;
889              else
890                preSolveFile=true;
891              break;
892            case PFI:
893              lpSolver->factorization()->setForrestTomlin(action==0);
894              break;
895            case CRASH:
896              doCrash=action;
897              break;
898            case MESSAGES:
899              lpSolver->messageHandler()->setPrefix(action!=0);
900              break;
901            case CHOLESKY:
902              choleskyType = action;
903              break;
904            case GAMMA:
905              gamma=action;
906              break;
907            case BARRIERSCALE:
908              scaleBarrier=action;
909              break;
910            case KKT:
911              doKKT=action;
912              break;
913            case CROSSOVER:
914              crossover=action;
915              break;
916            case GOMORYCUTS:
917              defaultSettings=false; // user knows what she is doing
918              gomoryAction = action;
919              break;
920            case PROBINGCUTS:
921              defaultSettings=false; // user knows what she is doing
922              probingAction = action;
923              break;
924            case KNAPSACKCUTS:
925              defaultSettings=false; // user knows what she is doing
926              knapsackAction = action;
927              break;
928            case REDSPLITCUTS:
929              defaultSettings=false; // user knows what she is doing
930              redsplitAction = action;
931              break;
932            case CLIQUECUTS:
933              defaultSettings=false; // user knows what she is doing
934              cliqueAction = action;
935              break;
936            case FLOWCUTS:
937              defaultSettings=false; // user knows what she is doing
938              flowAction = action;
939              break;
940            case MIXEDCUTS:
941              defaultSettings=false; // user knows what she is doing
942              mixedAction = action;
943              break;
944            case TWOMIRCUTS:
945              defaultSettings=false; // user knows what she is doing
946              twomirAction = action;
947              break;
948            case ROUNDING:
949              defaultSettings=false; // user knows what she is doing
950              useRounding = action;
951              break;
952            case FPUMP:
953              defaultSettings=false; // user knows what she is doing
954              useFpump=action;
955              break;
956            case CUTSSTRATEGY:
957              gomoryAction = action;
958              probingAction = action;
959              knapsackAction = action;
960              //redsplitAction = action;
961              cliqueAction = action;
962              flowAction = action;
963              mixedAction = action;
964              twomirAction = action;
965              parameters[whichParam(GOMORYCUTS,numberParameters,parameters)].setCurrentOption(action);
966              parameters[whichParam(PROBINGCUTS,numberParameters,parameters)].setCurrentOption(action);
967              parameters[whichParam(KNAPSACKCUTS,numberParameters,parameters)].setCurrentOption(action);
968              //parameters[whichParam(REDSPLITCUTS,numberParameters,parameters)].setCurrentOption(action);
969              parameters[whichParam(CLIQUECUTS,numberParameters,parameters)].setCurrentOption(action);
970              parameters[whichParam(FLOWCUTS,numberParameters,parameters)].setCurrentOption(action);
971              parameters[whichParam(MIXEDCUTS,numberParameters,parameters)].setCurrentOption(action);
972              parameters[whichParam(TWOMIRCUTS,numberParameters,parameters)].setCurrentOption(action);
973              break;
974            case HEURISTICSTRATEGY:
975              useRounding = action;
976              useGreedy = action;
977              useCombine = action;
978              //useLocalTree = action;
979              useFpump=action;
980              parameters[whichParam(ROUNDING,numberParameters,parameters)].setCurrentOption(action);
981              parameters[whichParam(GREEDY,numberParameters,parameters)].setCurrentOption(action);
982              parameters[whichParam(COMBINE,numberParameters,parameters)].setCurrentOption(action);
983              //parameters[whichParam(LOCALTREE,numberParameters,parameters)].setCurrentOption(action);
984              parameters[whichParam(FPUMP,numberParameters,parameters)].setCurrentOption(action);
985              break;
986            case GREEDY:
987              defaultSettings=false; // user knows what she is doing
988              useGreedy = action;
989              break;
990            case COMBINE:
991              defaultSettings=false; // user knows what she is doing
992              useCombine = action;
993              break;
994            case LOCALTREE:
995              defaultSettings=false; // user knows what she is doing
996              useLocalTree = action;
997              break;
998            case COSTSTRATEGY:
999              if (action!=1) {
1000                printf("Pseudo costs not implemented yet\n");
1001              } else {
1002                int numberColumns = model.getNumCols();
1003                int * sort = new int[numberColumns];
1004                double * dsort = new double[numberColumns];
1005                int * priority = new int [numberColumns];
1006                const double * objective = model.getObjCoefficients();
1007                int iColumn;
1008                int n=0;
1009                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1010                  if (model.isInteger(iColumn)) {
1011                    sort[n]=n;
1012                    dsort[n++]=-objective[iColumn];
1013                  }
1014                }
1015                CoinSort_2(dsort,dsort+n,sort);
1016                int level=0;
1017                double last = -1.0e100;
1018                for (int i=0;i<n;i++) {
1019                  int iPut=sort[i];
1020                  if (dsort[i]!=last) {
1021                    level++;
1022                    last=dsort[i];
1023                  }
1024                  priority[iPut]=level;
1025                }
1026                model.passInPriorities( priority,false);
1027                delete [] priority;
1028                delete [] sort;
1029                delete [] dsort;
1030              }
1031              break;
1032            case PREPROCESS:
1033              preProcess = action;
1034              break;
1035            default:
1036              abort();
1037            }
1038          }
1039        } else {
1040          // action
1041          if (type==EXIT)
1042            break; // stop all
1043          switch (type) {
1044          case DUALSIMPLEX:
1045          case PRIMALSIMPLEX:
1046          case SOLVECONTINUOUS:
1047          case BARRIER:
1048            if (goodModel) {
1049              ClpSolve::SolveType method;
1050              ClpSolve::PresolveType presolveType;
1051              ClpSimplex * model2 = lpSolver;
1052              ClpSolve solveOptions;
1053              if (preSolve!=5&&preSolve)
1054                presolveType=ClpSolve::presolveNumber;
1055              else if (preSolve)
1056                presolveType=ClpSolve::presolveOn;
1057              else
1058                presolveType=ClpSolve::presolveOff;
1059              solveOptions.setPresolveType(presolveType,preSolve);
1060              if (type==DUALSIMPLEX||SOLVECONTINUOUS) {
1061                method=ClpSolve::useDual;
1062              } else if (type==PRIMALSIMPLEX) {
1063                method=ClpSolve::usePrimalorSprint;
1064              } else {
1065                method = ClpSolve::useBarrier;
1066                if (crossover==1) {
1067                  method=ClpSolve::useBarrierNoCross;
1068                } else if (crossover==2) {
1069                  ClpObjective * obj = lpSolver->objectiveAsObject();
1070                  if (obj->type()>1) {
1071                    method=ClpSolve::useBarrierNoCross;
1072                    presolveType=ClpSolve::presolveOff;
1073                    solveOptions.setPresolveType(presolveType,preSolve);
1074                  } 
1075                }
1076              }
1077              solveOptions.setSolveType(method);
1078              if(preSolveFile)
1079                presolveOptions |= 0x40000000;
1080              solveOptions.setSpecialOption(4,presolveOptions);
1081              solveOptions.setSpecialOption(5,printOptions);
1082              if (method==ClpSolve::useDual) {
1083                // dual
1084                if (doCrash)
1085                  solveOptions.setSpecialOption(0,1,doCrash); // crash
1086                else if (doIdiot)
1087                  solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
1088              } else if (method==ClpSolve::usePrimalorSprint) {
1089                // primal
1090                // if slp turn everything off
1091                if (slpValue>0) {
1092                  doCrash=false;
1093                  doSprint=0;
1094                  doIdiot=-1;
1095                  solveOptions.setSpecialOption(1,10,slpValue); // slp
1096                  method=ClpSolve::usePrimal;
1097                }
1098                if (doCrash) {
1099                  solveOptions.setSpecialOption(1,1,doCrash); // crash
1100                } else if (doSprint>0) {
1101                  // sprint overrides idiot
1102                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
1103                } else if (doIdiot>0) {
1104                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
1105                } else if (slpValue<=0) {
1106                  if (doIdiot==0) {
1107                    if (doSprint==0)
1108                      solveOptions.setSpecialOption(1,4); // all slack
1109                    else
1110                      solveOptions.setSpecialOption(1,9); // all slack or sprint
1111                  } else {
1112                    if (doSprint==0)
1113                      solveOptions.setSpecialOption(1,8); // all slack or idiot
1114                    else
1115                      solveOptions.setSpecialOption(1,7); // initiative
1116                  }
1117                }
1118                if (basisHasValues==-1)
1119                  solveOptions.setSpecialOption(1,11); // switch off values
1120              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
1121                int barrierOptions = choleskyType;
1122                if (scaleBarrier)
1123                  barrierOptions |= 8;
1124                if (gamma)
1125                  barrierOptions |= 16;
1126                if (doKKT)
1127                  barrierOptions |= 32;
1128                solveOptions.setSpecialOption(4,barrierOptions);
1129              }
1130              model2->initialSolve(solveOptions);
1131              basisHasValues=1;
1132            } else {
1133              std::cout<<"** Current model not valid"<<std::endl;
1134            }
1135            break;
1136          case TIGHTEN:
1137            if (goodModel) {
1138              int numberInfeasibilities = lpSolver->tightenPrimalBounds();
1139              if (numberInfeasibilities)
1140                std::cout<<"** Analysis indicates model infeasible"<<std::endl;
1141            } else {
1142              std::cout<<"** Current model not valid"<<std::endl;
1143            }
1144            break;
1145          case PLUSMINUS:
1146            if (goodModel) {
1147              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
1148              ClpPackedMatrix* clpMatrix =
1149                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1150              if (clpMatrix) {
1151                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
1152                if (newMatrix->getIndices()) {
1153                  lpSolver->replaceMatrix(newMatrix);
1154                  delete saveMatrix;
1155                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
1156                } else {
1157                  std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
1158                }
1159              } else {
1160                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
1161              }
1162            } else {
1163              std::cout<<"** Current model not valid"<<std::endl;
1164            }
1165            break;
1166          case NETWORK:
1167            if (goodModel) {
1168              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
1169              ClpPackedMatrix* clpMatrix =
1170                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1171              if (clpMatrix) {
1172                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
1173                if (newMatrix->getIndices()) {
1174                  lpSolver->replaceMatrix(newMatrix);
1175                  delete saveMatrix;
1176                  std::cout<<"Matrix converted to network matrix"<<std::endl;
1177                } else {
1178                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
1179                }
1180              } else {
1181                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
1182              }
1183            } else {
1184              std::cout<<"** Current model not valid"<<std::endl;
1185            }
1186            break;
1187/*
1188  Run branch-and-cut. First set a few options -- node comparison, scaling. If
1189  the solver is Clp, consider running some presolve code (not yet converted
1190  this to generic OSI) with branch-and-cut. If presolve is disabled, or the
1191  solver is not Clp, simply run branch-and-cut. Print elapsed time at the end.
1192*/
1193          case BAB: // branchAndBound
1194          case STRENGTHEN:
1195            if (goodModel) {
1196              // If user made settings then use them
1197              if (!defaultSettings) {
1198                OsiSolverInterface * solver = model.solver();
1199                if (!doScaling)
1200                  solver->setHintParam(OsiDoScale,false,OsiHintTry);
1201                OsiClpSolverInterface * si =
1202                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
1203                assert (si != NULL);
1204                // get clp itself
1205                ClpSimplex * modelC = si->getModelPtr();
1206                if (modelC->tightenPrimalBounds()!=0) {
1207                  std::cout<<"Problem is infeasible!"<<std::endl;
1208                  break;
1209                }
1210                model.initialSolve();
1211                // bounds based on continuous
1212                if (tightenFactor) {
1213                  if (modelC->tightenPrimalBounds(tightenFactor)!=0) {
1214                    std::cout<<"Problem is infeasible!"<<std::endl;
1215                    break;
1216                  }
1217                }
1218                if (djFix<1.0e20) {
1219                  // do some fixing
1220                  int numberColumns = modelC->numberColumns();
1221                  int i;
1222                  const char * type = modelC->integerInformation();
1223                  double * lower = modelC->columnLower();
1224                  double * upper = modelC->columnUpper();
1225                  double * solution = modelC->primalColumnSolution();
1226                  double * dj = modelC->dualColumnSolution();
1227                  int numberFixed=0;
1228                  for (i=0;i<numberColumns;i++) {
1229                    if (type[i]) {
1230                      double value = solution[i];
1231                      if (value<lower[i]+1.0e-5&&dj[i]>djFix) {
1232                        solution[i]=lower[i];
1233                        upper[i]=lower[i];
1234                        numberFixed++;
1235                      } else if (value>upper[i]-1.0e-5&&dj[i]<-djFix) {
1236                        solution[i]=upper[i];
1237                        lower[i]=upper[i];
1238                        numberFixed++;
1239                      }
1240                    }
1241                  }
1242                  printf("%d columns fixed\n",numberFixed);
1243                }
1244              }
1245              // See if we want preprocessing
1246              OsiSolverInterface * saveSolver=NULL;
1247              CglPreProcess process;
1248              delete babModel;
1249              babModel = new CbcModel(model);
1250              OsiSolverInterface * solver3 = clpSolver->clone();
1251              if (defaultSettings) {
1252                OsiClpSolverInterface * si =
1253                  dynamic_cast<OsiClpSolverInterface *>(solver3) ;
1254                assert (si != NULL);
1255                // get clp itself
1256                ClpSimplex * modelC = si->getModelPtr();
1257                if (modelC->tightenPrimalBounds()!=0) {
1258                  std::cout<<"Problem is infeasible!"<<std::endl;
1259                  break;
1260                }
1261              }
1262              babModel->assignSolver(solver3);
1263              OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
1264              int numberChanged=0;
1265              if (clpSolver2->messageHandler()->logLevel())
1266                clpSolver2->messageHandler()->setLogLevel(1);
1267              int logLevel = parameters[slog].intValue();
1268              if (logLevel)
1269                clpSolver2->messageHandler()->setLogLevel(logLevel);
1270              lpSolver = clpSolver2->getModelPtr();
1271              if (lpSolver->factorizationFrequency()==200) {
1272                // User did not touch preset
1273                int numberRows = lpSolver->numberRows();
1274                const int cutoff1=10000;
1275                const int cutoff2=100000;
1276                const int base=75;
1277                const int freq0 = 50;
1278                const int freq1=200;
1279                const int freq2=400;
1280                const int maximum=1000;
1281                int frequency;
1282                if (numberRows<cutoff1)
1283                  frequency=base+numberRows/freq0;
1284                else if (numberRows<cutoff2)
1285                  frequency=base+cutoff1/freq0 + (numberRows-cutoff1)/freq1;
1286                else
1287                  frequency=base+cutoff1/freq0 + (cutoff2-cutoff1)/freq1 + (numberRows-cutoff2)/freq2;
1288                lpSolver->setFactorizationFrequency(CoinMin(maximum,frequency));
1289              }
1290              time2 = CoinCpuTime();
1291              totalTime += time2-time1;
1292              time1 = time2;
1293              double timeLeft = babModel->getMaximumSeconds();
1294              if (preProcess&&type==BAB) {
1295                saveSolver=babModel->solver()->clone();
1296                /* Do not try and produce equality cliques and
1297                   do up to 10 passes */
1298                OsiSolverInterface * solver2 = process.preProcess(*saveSolver,false,10);
1299                if (!solver2) {
1300                  printf("Pre-processing says infeasible\n");
1301                  break;
1302                } else {
1303                  printf("processed model has %d rows and %d columns\n",
1304                         solver2->getNumRows(),solver2->getNumCols());
1305                }
1306                //solver2->resolve();
1307                if (preProcess==2) {
1308                  OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (solver2);
1309                  ClpSimplex * lpSolver = clpSolver2->getModelPtr();
1310                  lpSolver->writeMps("presolved.mps",0,1,lpSolver->optimizationDirection());
1311                  printf("Preprocessed model (minimization) on presolved.mps\n");
1312                }
1313                // we have to keep solver2 so pass clone
1314                solver2 = solver2->clone();
1315                babModel->assignSolver(solver2);
1316                babModel->initialSolve();
1317                babModel->setMaximumSeconds(timeLeft-(CoinCpuTime()-time1));
1318              }
1319              if (0) {
1320                // for debug
1321                std::string problemName ;
1322                babModel->solver()->getStrParam(OsiProbName,problemName) ;
1323                babModel->solver()->activateRowCutDebugger(problemName.c_str()) ;
1324                twomirGen.probname_=strdup(problemName.c_str());
1325                // checking seems odd
1326                //redsplitGen.set_given_optsol(babModel->solver()->getRowCutDebuggerAlways()->optimalSolution(),
1327                //                         babModel->getNumCols());
1328              }
1329              if (debugValues) {
1330                if (numberDebugValues==babModel->getNumCols()) {
1331                  // for debug
1332                  babModel->solver()->activateRowCutDebugger(debugValues) ;
1333                } else {
1334                  printf("debug file has incorrect number of columns\n");
1335                }
1336              }
1337              // FPump done first as it only works if no solution
1338              CbcHeuristicFPump heuristic4(*babModel);
1339              if (useFpump) {
1340                heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
1341                babModel->addHeuristic(&heuristic4);
1342              }
1343              CbcRounding heuristic1(*babModel);
1344              if (useRounding)
1345                babModel->addHeuristic(&heuristic1) ;
1346              CbcHeuristicLocal heuristic2(*babModel);
1347              heuristic2.setSearchType(1);
1348              if (useCombine)
1349                babModel->addHeuristic(&heuristic2);
1350              CbcHeuristicGreedyCover heuristic3(*babModel);
1351              CbcHeuristicGreedyEquality heuristic3a(*babModel);
1352              if (useGreedy) {
1353                babModel->addHeuristic(&heuristic3);
1354                babModel->addHeuristic(&heuristic3a);
1355              }
1356              if (useLocalTree) {
1357                CbcTreeLocal localTree(babModel,NULL,10,0,0,10000,2000);
1358                babModel->passInTreeHandler(localTree);
1359              }
1360              // add cut generators if wanted
1361              if (probingAction==1)
1362                babModel->addCutGenerator(&probingGen,-1,"Probing");
1363              else if (probingAction>=2)
1364                babModel->addCutGenerator(&probingGen,-101+probingAction,"Probing");
1365              if (gomoryAction==1)
1366                babModel->addCutGenerator(&gomoryGen,-1,"Gomory");
1367              else if (gomoryAction>=2)
1368                babModel->addCutGenerator(&gomoryGen,-101+gomoryAction,"Gomory");
1369              if (knapsackAction==1)
1370                babModel->addCutGenerator(&knapsackGen,-1,"Knapsack");
1371              else if (knapsackAction>=2)
1372                babModel->addCutGenerator(&knapsackGen,-101+knapsackAction,"Knapsack");
1373              if (redsplitAction==1)
1374                babModel->addCutGenerator(&redsplitGen,-1,"Reduce-and-split");
1375              else if (redsplitAction>=2)
1376                babModel->addCutGenerator(&redsplitGen,-101+redsplitAction,"Reduce-and-split");
1377              if (cliqueAction==1)
1378                babModel->addCutGenerator(&cliqueGen,-1,"Clique");
1379              else if (cliqueAction>=2)
1380                babModel->addCutGenerator(&cliqueGen,-101+cliqueAction,"Clique");
1381              if (mixedAction==1)
1382                babModel->addCutGenerator(&mixedGen,-1,"MixedIntegerRounding2");
1383              else if (mixedAction>=2)
1384                babModel->addCutGenerator(&mixedGen,-101+mixedAction,"MixedIntegerRounding2");
1385              if (flowAction==1)
1386                babModel->addCutGenerator(&flowGen,-1,"FlowCover");
1387              else if (flowAction>=2)
1388                babModel->addCutGenerator(&flowGen,-101+flowAction,"FlowCover");
1389              if (twomirAction==1)
1390                babModel->addCutGenerator(&twomirGen,-1,"TwoMirCuts");
1391              else if (twomirAction>=2)
1392                babModel->addCutGenerator(&twomirGen,-101+twomirAction,"TwoMirCuts");
1393              // Say we want timings
1394              int numberGenerators = babModel->numberCutGenerators();
1395              int iGenerator;
1396              int cutDepth=
1397                parameters[whichParam(CUTDEPTH,numberParameters,parameters)].intValue();
1398              for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
1399                CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
1400                if (generator->howOften()<=-98)
1401                  generator->setSwitchOffIfLessThan(1);
1402                generator->setTiming(true);
1403                if (cutDepth>=0)
1404                  generator->setWhatDepth(cutDepth) ;
1405              }
1406              // Could tune more
1407              babModel->setMinimumDrop(min(5.0e-2,
1408                                        fabs(babModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
1409              if (cutPass==-1234567) {
1410                if (babModel->getNumCols()<500)
1411                  babModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
1412                else if (babModel->getNumCols()<5000)
1413                  babModel->setMaximumCutPassesAtRoot(100); // use minimum drop
1414                else
1415                  babModel->setMaximumCutPassesAtRoot(20);
1416              } else {
1417                babModel->setMaximumCutPassesAtRoot(cutPass);
1418              }
1419              babModel->setMaximumCutPasses(1);
1420             
1421              // Do more strong branching if small
1422              //if (babModel->getNumCols()<5000)
1423              //babModel->setNumberStrong(20);
1424              // Switch off strong branching if wanted
1425              //if (babModel->getNumCols()>10*babModel->getNumRows())
1426              //babModel->setNumberStrong(0);
1427              babModel->messageHandler()->setLogLevel(parameters[log].intValue());
1428              if (babModel->getNumCols()>2000||babModel->getNumRows()>1500||
1429                  babModel->messageHandler()->logLevel()>1)
1430                babModel->setPrintFrequency(100);
1431             
1432              babModel->solver()->setIntParam(OsiMaxNumIterationHotStart,100);
1433              OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
1434              // go faster stripes
1435              if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
1436                osiclp->setupForRepeatedUse(2,0);
1437              } else {
1438                osiclp->setupForRepeatedUse(0,0);
1439              }
1440              double increment=babModel->getCutoffIncrement();;
1441              int * changed = analyze( osiclp,numberChanged,increment,false);
1442              babModel->setCutoffIncrement(CoinMax(babModel->getCutoffIncrement(),increment));
1443              // Turn this off if you get problems
1444              // Used to be automatically set
1445              osiclp->setSpecialOptions(osiclp->specialOptions()|(128+64));
1446              if (gapRatio < 1.0e100) {
1447                double value = babModel->solver()->getObjValue() ;
1448                double value2 = gapRatio*(1.0e-5+fabs(value)) ;
1449                babModel->setAllowableGap(value2) ;
1450                std::cout << "Continuous " << value
1451                          << ", so allowable gap set to "
1452                          << value2 << std::endl ;
1453              }
1454              // probably faster to use a basis to get integer solutions
1455              babModel->setSpecialOptions(2);
1456              currentBranchModel = babModel;
1457              OsiSolverInterface * strengthenedModel=NULL;
1458              if (type==BAB) {
1459                babModel->branchAndBound();
1460              } else {
1461                strengthenedModel = babModel->strengthenedModel();
1462              }
1463              currentBranchModel = NULL;
1464              osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
1465              if (debugFile=="createAfterPre"&&babModel->bestSolution()) {
1466                lpSolver = osiclp->getModelPtr();
1467                //move best solution (should be there -- but ..)
1468                int n = lpSolver->getNumCols();
1469                memcpy(lpSolver->primalColumnSolution(),babModel->bestSolution(),n*sizeof(double));
1470                saveSolution(osiclp->getModelPtr(),"debug.file");
1471              }
1472              // Print more statistics
1473              std::cout<<"Cuts at root node changed objective from "<<babModel->getContinuousObjective()
1474                       <<" to "<<babModel->rootObjectiveAfterCuts()<<std::endl;
1475             
1476              for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
1477                CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
1478                std::cout<<generator->cutGeneratorName()<<" was tried "
1479                         <<generator->numberTimesEntered()<<" times and created "
1480                         <<generator->numberCutsInTotal()<<" cuts of which "
1481                         <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
1482                if (generator->timing())
1483                  std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
1484                else
1485                  std::cout<<std::endl;
1486              }
1487              time2 = CoinCpuTime();
1488              totalTime += time2-time1;
1489              if (babModel->getMinimizationObjValue()<1.0e50&&type==BAB) {
1490                // post process
1491                if (preProcess) {
1492                  process.postProcess(*babModel->solver());
1493                  // Solution now back in saveSolver
1494                  babModel->assignSolver(saveSolver);
1495                }
1496              }
1497              if (type==STRENGTHEN&&strengthenedModel)
1498                clpSolver = dynamic_cast< OsiClpSolverInterface*> (strengthenedModel);
1499              lpSolver = clpSolver->getModelPtr();
1500              if (debugFile=="create"&&babModel->bestSolution()) {
1501                saveSolution(lpSolver,"debug.file");
1502              }
1503              if (numberChanged) {
1504                for (int i=0;i<numberChanged;i++) {
1505                  int iColumn=changed[i];
1506                  clpSolver->setContinuous(iColumn);
1507                }
1508                delete [] changed;
1509              }
1510              if (type==BAB) {
1511                std::string statusName[]={"Finished","Stopped on ","Difficulties",
1512                                          "","","User ctrl-c"};
1513                std::string minor[]={"","","gap","nodes","time","","solutions"};
1514                int iStat = babModel->status();
1515                int iStat2 = babModel->secondaryStatus();
1516                std::cout<<"Result - "<<statusName[iStat]<<minor[iStat2]
1517                         <<" objective "<<babModel->getObjValue()<<
1518                  " after "<<babModel->getNodeCount()<<" nodes and "
1519                         <<babModel->getIterationCount()<<
1520                  " iterations - took "<<time2-time1<<" seconds"<<std::endl;
1521              } else {
1522                std::cout<<"Model strengthend - now has "<<clpSolver->getNumRows()
1523                         <<" rows"<<std::endl;
1524              }
1525              time1 = time2;
1526              delete babModel;
1527              babModel=NULL;
1528            } else {
1529              std::cout << "** Current model not valid" << std::endl ; 
1530            }
1531            break ;
1532          case IMPORT:
1533            {
1534              delete babModel;
1535              babModel=NULL;
1536              // get next field
1537              field = CoinReadGetString(argc,argv);
1538              if (field=="$") {
1539                field = parameters[iParam].stringValue();
1540              } else if (field=="EOL") {
1541                parameters[iParam].printString();
1542                break;
1543              } else {
1544                parameters[iParam].setStringValue(field);
1545              }
1546              std::string fileName;
1547              bool canOpen=false;
1548              if (field=="-") {
1549                // stdin
1550                canOpen=true;
1551                fileName = "-";
1552              } else {
1553                bool absolutePath;
1554                if (dirsep=='/') {
1555                  // non Windows (or cygwin)
1556                  absolutePath=(field[0]=='/');
1557                } else {
1558                  //Windows (non cycgwin)
1559                  absolutePath=(field[0]=='\\');
1560                  // but allow for :
1561                  if (strchr(field.c_str(),':'))
1562                    absolutePath=true;
1563                }
1564                if (absolutePath) {
1565                  fileName = field;
1566                } else if (field[0]=='~') {
1567                  char * environ = getenv("HOME");
1568                  if (environ) {
1569                    std::string home(environ);
1570                    field=field.erase(0,1);
1571                    fileName = home+field;
1572                  } else {
1573                    fileName=field;
1574                  }
1575                } else {
1576                  fileName = directory+field;
1577                }
1578                FILE *fp=fopen(fileName.c_str(),"r");
1579                if (fp) {
1580                  // can open - lets go for it
1581                  fclose(fp);
1582                  canOpen=true;
1583                } else {
1584                  std::cout<<"Unable to open file "<<fileName<<std::endl;
1585                }
1586              }
1587              if (canOpen) {
1588                int status =lpSolver->readMps(fileName.c_str(),
1589                                                   keepImportNames!=0,
1590                                                   allowImportErrors!=0);
1591                if (!status||(status>0&&allowImportErrors)) {
1592                  if (keepImportNames) {
1593                    lengthName = lpSolver->lengthNames();
1594                    rowNames = *(lpSolver->rowNames());
1595                    columnNames = *(lpSolver->columnNames());
1596                  } else {
1597                    lengthName=0;
1598                  }
1599                  goodModel=true;
1600                  //Set integers in clpsolver
1601                  const char * info = lpSolver->integerInformation();
1602                  if (info) {
1603                    int numberColumns = lpSolver->numberColumns();
1604                    int i;
1605                    for (i=0;i<numberColumns;i++) {
1606                      if (info[i]) 
1607                        clpSolver->setInteger(i);
1608                    }
1609                  }
1610                  // sets to all slack (not necessary?)
1611                  lpSolver->createStatus();
1612                  time2 = CoinCpuTime();
1613                  totalTime += time2-time1;
1614                  time1=time2;
1615                  // Go to canned file if just input file
1616                  if (CbcOrClpRead_mode==2&&argc==2) {
1617                    // only if ends .mps
1618                    char * find = strstr(fileName.c_str(),".mps");
1619                    if (find&&find[4]=='\0') {
1620                      find[1]='p'; find[2]='a';find[3]='r';
1621                      FILE *fp=fopen(fileName.c_str(),"r");
1622                      if (fp) {
1623                        CbcOrClpReadCommand=fp; // Read from that file
1624                        CbcOrClpRead_mode=-1;
1625                      }
1626                    }
1627                  }
1628                } else {
1629                  // errors
1630                  std::cout<<"There were "<<status<<
1631                    " errors on input"<<std::endl;
1632                }
1633              }
1634            }
1635            break;
1636          case EXPORT:
1637            if (goodModel) {
1638              // get next field
1639              field = CoinReadGetString(argc,argv);
1640              if (field=="$") {
1641                field = parameters[iParam].stringValue();
1642              } else if (field=="EOL") {
1643                parameters[iParam].printString();
1644                break;
1645              } else {
1646                parameters[iParam].setStringValue(field);
1647              }
1648              std::string fileName;
1649              bool canOpen=false;
1650              if (field[0]=='/'||field[0]=='\\') {
1651                fileName = field;
1652              } else if (field[0]=='~') {
1653                char * environ = getenv("HOME");
1654                if (environ) {
1655                  std::string home(environ);
1656                  field=field.erase(0,1);
1657                  fileName = home+field;
1658                } else {
1659                  fileName=field;
1660                }
1661              } else {
1662                fileName = directory+field;
1663              }
1664              FILE *fp=fopen(fileName.c_str(),"w");
1665              if (fp) {
1666                // can open - lets go for it
1667                fclose(fp);
1668                canOpen=true;
1669              } else {
1670                std::cout<<"Unable to open file "<<fileName<<std::endl;
1671              }
1672              if (canOpen) {
1673                // If presolve on then save presolved
1674                bool deleteModel2=false;
1675                ClpSimplex * model2 = lpSolver;
1676                if (preSolve) {
1677                  ClpPresolve pinfo;
1678                  int presolveOptions2 = presolveOptions&~0x40000000;
1679                  if ((presolveOptions2&0xffff)!=0)
1680                    pinfo.setPresolveActions(presolveOptions2);
1681                  if ((printOptions&1)!=0)
1682                    pinfo.statistics();
1683                  model2 = 
1684                    pinfo.presolvedModel(*lpSolver,1.0e-8,
1685                                         true,preSolve);
1686                  if (model2) {
1687                    printf("Saving presolved model on %s\n",
1688                           fileName.c_str());
1689                    deleteModel2=true;
1690                  } else {
1691                    printf("Presolved model looks infeasible - saving original on %s\n",
1692                           fileName.c_str());
1693                    deleteModel2=false;
1694                    model2 = lpSolver;
1695
1696                  }
1697                } else {
1698                  printf("Saving model on %s\n",
1699                           fileName.c_str());
1700                }
1701#if 0
1702                // Convert names
1703                int iRow;
1704                int numberRows=model2->numberRows();
1705                int iColumn;
1706                int numberColumns=model2->numberColumns();
1707
1708                char ** rowNames = NULL;
1709                char ** columnNames = NULL;
1710                if (model2->lengthNames()) {
1711                  rowNames = new char * [numberRows];
1712                  for (iRow=0;iRow<numberRows;iRow++) {
1713                    rowNames[iRow] =
1714                      strdup(model2->rowName(iRow).c_str());
1715#ifdef STRIPBLANKS
1716                    char * xx = rowNames[iRow];
1717                    int i;
1718                    int length = strlen(xx);
1719                    int n=0;
1720                    for (i=0;i<length;i++) {
1721                      if (xx[i]!=' ')
1722                        xx[n++]=xx[i];
1723                    }
1724                    xx[n]='\0';
1725#endif
1726                  }
1727                 
1728                  columnNames = new char * [numberColumns];
1729                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1730                    columnNames[iColumn] =
1731                      strdup(model2->columnName(iColumn).c_str());
1732#ifdef STRIPBLANKS
1733                    char * xx = columnNames[iColumn];
1734                    int i;
1735                    int length = strlen(xx);
1736                    int n=0;
1737                    for (i=0;i<length;i++) {
1738                      if (xx[i]!=' ')
1739                        xx[n++]=xx[i];
1740                    }
1741                    xx[n]='\0';
1742#endif
1743                  }
1744                }
1745                CoinMpsIO writer;
1746                writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
1747                                  model2->getColLower(), model2->getColUpper(),
1748                                  model2->getObjCoefficients(),
1749                                  (const char*) 0 /*integrality*/,
1750                                  model2->getRowLower(), model2->getRowUpper(),
1751                                  columnNames, rowNames);
1752                // Pass in array saying if each variable integer
1753                writer.copyInIntegerInformation(model2->integerInformation());
1754                writer.setObjectiveOffset(model2->objectiveOffset());
1755                writer.writeMps(fileName.c_str(),0,1,1);
1756                if (rowNames) {
1757                  for (iRow=0;iRow<numberRows;iRow++) {
1758                    free(rowNames[iRow]);
1759                  }
1760                  delete [] rowNames;
1761                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1762                    free(columnNames[iColumn]);
1763                  }
1764                  delete [] columnNames;
1765                }
1766#else
1767                model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
1768#endif
1769                if (deleteModel2)
1770                  delete model2;
1771                time2 = CoinCpuTime();
1772                totalTime += time2-time1;
1773                time1=time2;
1774              }
1775            } else {
1776              std::cout<<"** Current model not valid"<<std::endl;
1777            }
1778            break;
1779          case BASISIN:
1780            if (goodModel) {
1781              // get next field
1782              field = CoinReadGetString(argc,argv);
1783              if (field=="$") {
1784                field = parameters[iParam].stringValue();
1785              } else if (field=="EOL") {
1786                parameters[iParam].printString();
1787                break;
1788              } else {
1789                parameters[iParam].setStringValue(field);
1790              }
1791              std::string fileName;
1792              bool canOpen=false;
1793              if (field=="-") {
1794                // stdin
1795                canOpen=true;
1796                fileName = "-";
1797              } else {
1798                if (field[0]=='/'||field[0]=='\\') {
1799                  fileName = field;
1800                } else if (field[0]=='~') {
1801                  char * environ = getenv("HOME");
1802                  if (environ) {
1803                    std::string home(environ);
1804                    field=field.erase(0,1);
1805                    fileName = home+field;
1806                  } else {
1807                    fileName=field;
1808                  }
1809                } else {
1810                  fileName = directory+field;
1811                }
1812                FILE *fp=fopen(fileName.c_str(),"r");
1813                if (fp) {
1814                  // can open - lets go for it
1815                  fclose(fp);
1816                  canOpen=true;
1817                } else {
1818                  std::cout<<"Unable to open file "<<fileName<<std::endl;
1819                }
1820              }
1821              if (canOpen) {
1822                int values = lpSolver->readBasis(fileName.c_str());
1823                if (values==0)
1824                  basisHasValues=-1;
1825                else
1826                  basisHasValues=1;
1827              }
1828            } else {
1829              std::cout<<"** Current model not valid"<<std::endl;
1830            }
1831            break;
1832          case DEBUG:
1833            if (goodModel) {
1834              delete [] debugValues;
1835              debugValues=NULL;
1836              // get next field
1837              field = CoinReadGetString(argc,argv);
1838              if (field=="$") {
1839                field = parameters[iParam].stringValue();
1840              } else if (field=="EOL") {
1841                parameters[iParam].printString();
1842                break;
1843              } else {
1844                parameters[iParam].setStringValue(field);
1845                debugFile=field;
1846                if (debugFile=="create"||
1847                    debugFile=="createAfterPre") {
1848                  printf("Will create a debug file so this run should be a good one\n");
1849                  break;
1850                }
1851              }
1852              std::string fileName;
1853              if (field[0]=='/'||field[0]=='\\') {
1854                fileName = field;
1855              } else if (field[0]=='~') {
1856                char * environ = getenv("HOME");
1857                if (environ) {
1858                  std::string home(environ);
1859                  field=field.erase(0,1);
1860                  fileName = home+field;
1861                } else {
1862                  fileName=field;
1863                }
1864              } else {
1865                fileName = directory+field;
1866              }
1867              FILE *fp=fopen(fileName.c_str(),"rb");
1868              if (fp) {
1869                // can open - lets go for it
1870                int numRows;
1871                double obj;
1872                fread(&numRows,sizeof(int),1,fp);
1873                fread(&numberDebugValues,sizeof(int),1,fp);
1874                fread(&obj,sizeof(double),1,fp);
1875                debugValues = new double[numberDebugValues+numRows];
1876                fread(debugValues,sizeof(double),numRows,fp);
1877                fread(debugValues,sizeof(double),numRows,fp);
1878                fread(debugValues,sizeof(double),numberDebugValues,fp);
1879                printf("%d doubles read into debugValues\n",numberDebugValues);
1880                fclose(fp);
1881              } else {
1882                std::cout<<"Unable to open file "<<fileName<<std::endl;
1883              }
1884            } else {
1885              std::cout<<"** Current model not valid"<<std::endl;
1886            }
1887            break;
1888          case BASISOUT:
1889            if (goodModel) {
1890              // get next field
1891              field = CoinReadGetString(argc,argv);
1892              if (field=="$") {
1893                field = parameters[iParam].stringValue();
1894              } else if (field=="EOL") {
1895                parameters[iParam].printString();
1896                break;
1897              } else {
1898                parameters[iParam].setStringValue(field);
1899              }
1900              std::string fileName;
1901              bool canOpen=false;
1902              if (field[0]=='/'||field[0]=='\\') {
1903                fileName = field;
1904              } else if (field[0]=='~') {
1905                char * environ = getenv("HOME");
1906                if (environ) {
1907                  std::string home(environ);
1908                  field=field.erase(0,1);
1909                  fileName = home+field;
1910                } else {
1911                  fileName=field;
1912                }
1913              } else {
1914                fileName = directory+field;
1915              }
1916              FILE *fp=fopen(fileName.c_str(),"w");
1917              if (fp) {
1918                // can open - lets go for it
1919                fclose(fp);
1920                canOpen=true;
1921              } else {
1922                std::cout<<"Unable to open file "<<fileName<<std::endl;
1923              }
1924              if (canOpen) {
1925                ClpSimplex * model2 = lpSolver;
1926                model2->writeBasis(fileName.c_str(),outputFormat>1,outputFormat-2);
1927                time2 = CoinCpuTime();
1928                totalTime += time2-time1;
1929                time1=time2;
1930              }
1931            } else {
1932              std::cout<<"** Current model not valid"<<std::endl;
1933            }
1934            break;
1935          case SAVE:
1936            {
1937              // get next field
1938              field = CoinReadGetString(argc,argv);
1939              if (field=="$") {
1940                field = parameters[iParam].stringValue();
1941              } else if (field=="EOL") {
1942                parameters[iParam].printString();
1943                break;
1944              } else {
1945                parameters[iParam].setStringValue(field);
1946              }
1947              std::string fileName;
1948              bool canOpen=false;
1949              if (field[0]=='/'||field[0]=='\\') {
1950                fileName = field;
1951              } else if (field[0]=='~') {
1952                char * environ = getenv("HOME");
1953                if (environ) {
1954                  std::string home(environ);
1955                  field=field.erase(0,1);
1956                  fileName = home+field;
1957                } else {
1958                  fileName=field;
1959                }
1960              } else {
1961                fileName = directory+field;
1962              }
1963              FILE *fp=fopen(fileName.c_str(),"wb");
1964              if (fp) {
1965                // can open - lets go for it
1966                fclose(fp);
1967                canOpen=true;
1968              } else {
1969                std::cout<<"Unable to open file "<<fileName<<std::endl;
1970              }
1971              if (canOpen) {
1972                int status;
1973                // If presolve on then save presolved
1974                bool deleteModel2=false;
1975                ClpSimplex * model2 = lpSolver;
1976                if (preSolve) {
1977                  ClpPresolve pinfo;
1978                  model2 = 
1979                    pinfo.presolvedModel(*lpSolver,1.0e-8,
1980                                         false,preSolve);
1981                  if (model2) {
1982                    printf("Saving presolved model on %s\n",
1983                           fileName.c_str());
1984                    deleteModel2=true;
1985                  } else {
1986                    printf("Presolved model looks infeasible - saving original on %s\n",
1987                           fileName.c_str());
1988                    deleteModel2=false;
1989                    model2 = lpSolver;
1990
1991                  }
1992                } else {
1993                  printf("Saving model on %s\n",
1994                           fileName.c_str());
1995                }
1996                status =model2->saveModel(fileName.c_str());
1997                if (deleteModel2)
1998                  delete model2;
1999                if (!status) {
2000                  goodModel=true;
2001                  time2 = CoinCpuTime();
2002                  totalTime += time2-time1;
2003                  time1=time2;
2004                } else {
2005                  // errors
2006                  std::cout<<"There were errors on output"<<std::endl;
2007                }
2008              }
2009            }
2010            break;
2011          case RESTORE:
2012            {
2013              // get next field
2014              field = CoinReadGetString(argc,argv);
2015              if (field=="$") {
2016                field = parameters[iParam].stringValue();
2017              } else if (field=="EOL") {
2018                parameters[iParam].printString();
2019                break;
2020              } else {
2021                parameters[iParam].setStringValue(field);
2022              }
2023              std::string fileName;
2024              bool canOpen=false;
2025              if (field[0]=='/'||field[0]=='\\') {
2026                fileName = field;
2027              } else if (field[0]=='~') {
2028                char * environ = getenv("HOME");
2029                if (environ) {
2030                  std::string home(environ);
2031                  field=field.erase(0,1);
2032                  fileName = home+field;
2033                } else {
2034                  fileName=field;
2035                }
2036              } else {
2037                fileName = directory+field;
2038              }
2039              FILE *fp=fopen(fileName.c_str(),"rb");
2040              if (fp) {
2041                // can open - lets go for it
2042                fclose(fp);
2043                canOpen=true;
2044              } else {
2045                std::cout<<"Unable to open file "<<fileName<<std::endl;
2046              }
2047              if (canOpen) {
2048                int status =lpSolver->restoreModel(fileName.c_str());
2049                if (!status) {
2050                  goodModel=true;
2051                  time2 = CoinCpuTime();
2052                  totalTime += time2-time1;
2053                  time1=time2;
2054                } else {
2055                  // errors
2056                  std::cout<<"There were errors on input"<<std::endl;
2057                }
2058              }
2059            }
2060            break;
2061          case MAXIMIZE:
2062            lpSolver->setOptimizationDirection(-1);
2063            break;
2064          case MINIMIZE:
2065            lpSolver->setOptimizationDirection(1);
2066            break;
2067          case ALLSLACK:
2068            lpSolver->allSlackBasis(true);
2069            break;
2070          case REVERSE:
2071            if (goodModel) {
2072              int iColumn;
2073              int numberColumns=lpSolver->numberColumns();
2074              double * dualColumnSolution = 
2075                lpSolver->dualColumnSolution();
2076              ClpObjective * obj = lpSolver->objectiveAsObject();
2077              assert(dynamic_cast<ClpLinearObjective *> (obj));
2078              double offset;
2079              double * objective = obj->gradient(NULL,NULL,offset,true);
2080              for (iColumn=0;iColumn<numberColumns;iColumn++) {
2081                dualColumnSolution[iColumn] = dualColumnSolution[iColumn];
2082                objective[iColumn] = -objective[iColumn];
2083              }
2084              int iRow;
2085              int numberRows=lpSolver->numberRows();
2086              double * dualRowSolution = 
2087                lpSolver->dualRowSolution();
2088              for (iRow=0;iRow<numberRows;iRow++) 
2089                dualRowSolution[iRow] = dualRowSolution[iRow];
2090            }
2091            break;
2092          case DIRECTORY:
2093            {
2094              std::string name = CoinReadGetString(argc,argv);
2095              if (name!="EOL") {
2096                int length=name.length();
2097                if (name[length-1]=='/'||name[length-1]=='\\')
2098                  directory=name;
2099                else
2100                  directory = name+"/";
2101                parameters[iParam].setStringValue(directory);
2102              } else {
2103                parameters[iParam].printString();
2104              }
2105            }
2106            break;
2107          case STDIN:
2108            CbcOrClpRead_mode=-1;
2109            break;
2110          case NETLIB_DUAL:
2111          case NETLIB_EITHER:
2112          case NETLIB_BARRIER:
2113          case NETLIB_PRIMAL:
2114          case NETLIB_TUNE:
2115            {
2116              // create fields for unitTest
2117              const char * fields[4];
2118              int nFields=2;
2119              fields[0]="fake main from unitTest";
2120              fields[1]="-netlib";
2121              if (directory!=defaultDirectory) {
2122                fields[2]="-netlibDir";
2123                fields[3]=directory.c_str();
2124                nFields=4;
2125              }
2126              int algorithm;
2127              if (type==NETLIB_DUAL) {
2128                std::cerr<<"Doing netlib with dual algorithm"<<std::endl;
2129                algorithm =0;
2130              } else if (type==NETLIB_BARRIER) {
2131                std::cerr<<"Doing netlib with barrier algorithm"<<std::endl;
2132                algorithm =2;
2133              } else if (type==NETLIB_EITHER) {
2134                std::cerr<<"Doing netlib with dual or primal algorithm"<<std::endl;
2135                algorithm =3;
2136              } else if (type==NETLIB_TUNE) {
2137                std::cerr<<"Doing netlib with best algorithm!"<<std::endl;
2138                algorithm =5;
2139                // uncomment next to get active tuning
2140                // algorithm=6;
2141              } else {
2142                std::cerr<<"Doing netlib with primal agorithm"<<std::endl;
2143                algorithm=1;
2144              }
2145              int specialOptions = lpSolver->specialOptions();
2146              lpSolver->setSpecialOptions(0);
2147              mainTest(nFields,fields,algorithm,*lpSolver,
2148                       (preSolve!=0),specialOptions);
2149            }
2150            break;
2151          case UNITTEST:
2152            {
2153              // create fields for unitTest
2154              const char * fields[3];
2155              int nFields=1;
2156              fields[0]="fake main from unitTest";
2157              if (directory!=defaultDirectory) {
2158                fields[1]="-mpsDir";
2159                fields[2]=directory.c_str();
2160                nFields=3;
2161              }
2162              mainTest(nFields,fields,false,*lpSolver,(preSolve!=0),
2163                       false);
2164            }
2165            break;
2166          case FAKEBOUND:
2167            if (goodModel) {
2168              // get bound
2169              double value = CoinReadGetDoubleField(argc,argv,&valid);
2170              if (!valid) {
2171                std::cout<<"Setting "<<parameters[iParam].name()<<
2172                  " to DEBUG "<<value<<std::endl;
2173                int iRow;
2174                int numberRows=lpSolver->numberRows();
2175                double * rowLower = lpSolver->rowLower();
2176                double * rowUpper = lpSolver->rowUpper();
2177                for (iRow=0;iRow<numberRows;iRow++) {
2178                  // leave free ones for now
2179                  if (rowLower[iRow]>-1.0e20||rowUpper[iRow]<1.0e20) {
2180                    rowLower[iRow]=CoinMax(rowLower[iRow],-value);
2181                    rowUpper[iRow]=CoinMin(rowUpper[iRow],value);
2182                  }
2183                }
2184                int iColumn;
2185                int numberColumns=lpSolver->numberColumns();
2186                double * columnLower = lpSolver->columnLower();
2187                double * columnUpper = lpSolver->columnUpper();
2188                for (iColumn=0;iColumn<numberColumns;iColumn++) {
2189                  // leave free ones for now
2190                  if (columnLower[iColumn]>-1.0e20||
2191                      columnUpper[iColumn]<1.0e20) {
2192                    columnLower[iColumn]=CoinMax(columnLower[iColumn],-value);
2193                    columnUpper[iColumn]=CoinMin(columnUpper[iColumn],value);
2194                  }
2195                }
2196              } else if (valid==1) {
2197                abort();
2198              } else {
2199                std::cout<<"enter value for "<<parameters[iParam].name()<<
2200                  std::endl;
2201              }
2202            }
2203            break;
2204          case REALLY_SCALE:
2205            if (goodModel) {
2206              ClpSimplex newModel(*lpSolver,
2207                                  lpSolver->scalingFlag());
2208              printf("model really really scaled\n");
2209              *lpSolver=newModel;
2210            }
2211            break;
2212          case HELP:
2213            std::cout<<"Coin Solver version "<<CBCVERSION
2214                     <<", build "<<__DATE__<<std::endl;
2215            std::cout<<"Non default values:-"<<std::endl;
2216            std::cout<<"Perturbation "<<lpSolver->perturbation()<<" (default 100)"
2217                     <<std::endl;
2218            CoinReadPrintit(
2219                    "Presolve being done with 5 passes\n\
2220Dual steepest edge steep/partial on matrix shape and factorization density\n\
2221Clpnnnn taken out of messages\n\
2222If Factorization frequency default then done on size of matrix\n\n\
2223(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
2224You can switch to interactive mode at any time so\n\
2225clp watson.mps -scaling off -primalsimplex\nis the same as\n\
2226clp watson.mps -\nscaling off\nprimalsimplex"
2227                    );
2228            break;
2229          case SOLUTION:
2230            if (goodModel) {
2231              // get next field
2232              field = CoinReadGetString(argc,argv);
2233              if (field=="$") {
2234                field = parameters[iParam].stringValue();
2235              } else if (field=="EOL") {
2236                parameters[iParam].printString();
2237                break;
2238              } else {
2239                parameters[iParam].setStringValue(field);
2240              }
2241              std::string fileName;
2242              FILE *fp=NULL;
2243              if (field=="-"||field=="EOL"||field=="stdout") {
2244                // stdout
2245                fp=stdout;
2246              } else if (field=="stderr") {
2247                // stderr
2248                fp=stderr;
2249              } else {
2250                if (field[0]=='/'||field[0]=='\\') {
2251                  fileName = field;
2252                } else if (field[0]=='~') {
2253                  char * environ = getenv("HOME");
2254                  if (environ) {
2255                    std::string home(environ);
2256                    field=field.erase(0,1);
2257                    fileName = home+field;
2258                  } else {
2259                    fileName=field;
2260                  }
2261                } else {
2262                  fileName = directory+field;
2263                }
2264                fp=fopen(fileName.c_str(),"w");
2265              }
2266              if (fp) {
2267                // make fancy later on
2268                int iRow;
2269                int numberRows=lpSolver->numberRows();
2270                double * dualRowSolution = lpSolver->dualRowSolution();
2271                double * primalRowSolution = 
2272                  lpSolver->primalRowSolution();
2273                double * rowLower = lpSolver->rowLower();
2274                double * rowUpper = lpSolver->rowUpper();
2275                double primalTolerance = lpSolver->primalTolerance();
2276                char format[6];
2277                sprintf(format,"%%-%ds",CoinMax(lengthName,8));
2278                if (printMode>2) {
2279                  for (iRow=0;iRow<numberRows;iRow++) {
2280                    int type=printMode-3;
2281                    if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
2282                        primalRowSolution[iRow]<rowLower[iRow]-primalTolerance) {
2283                      fprintf(fp,"** ");
2284                      type=2;
2285                    } else if (fabs(primalRowSolution[iRow])>1.0e-8) {
2286                      type=1;
2287                    } else if (numberRows<50) {
2288                      type=3;
2289                    }
2290                    if (type) {
2291                      fprintf(fp,"%7d ",iRow);
2292                      if (lengthName)
2293                        fprintf(fp,format,rowNames[iRow].c_str());
2294                      fprintf(fp,"%15.8g        %15.8g\n",primalRowSolution[iRow],
2295                              dualRowSolution[iRow]);
2296                    }
2297                  }
2298                }
2299                int iColumn;
2300                int numberColumns=lpSolver->numberColumns();
2301                double * dualColumnSolution = 
2302                  lpSolver->dualColumnSolution();
2303                double * primalColumnSolution = 
2304                  lpSolver->primalColumnSolution();
2305                double * columnLower = lpSolver->columnLower();
2306                double * columnUpper = lpSolver->columnUpper();
2307                if (printMode!=2) {
2308                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2309                    int type=0;
2310                    if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
2311                        primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance) {
2312                      fprintf(fp,"** ");
2313                      type=2;
2314                    } else if (fabs(primalColumnSolution[iColumn])>1.0e-8) {
2315                      type=1;
2316                    } else if (numberColumns<50) {
2317                      type=3;
2318                    }
2319                    // see if integer
2320                    if (!lpSolver->isInteger(iColumn)&&printMode==1)
2321                      type=0;
2322                    if (type) {
2323                      fprintf(fp,"%7d ",iColumn);
2324                      if (lengthName)
2325                        fprintf(fp,format,columnNames[iColumn].c_str());
2326                      fprintf(fp,"%15.8g        %15.8g\n",
2327                              primalColumnSolution[iColumn],
2328                              dualColumnSolution[iColumn]);
2329                    }
2330                  }
2331                } else {
2332                  // special format suitable for OsiRowCutDebugger
2333                  int n=0;
2334                  bool comma=false;
2335                  bool newLine=false;
2336                  fprintf(fp,"\tint intIndicesV[]={\n");
2337                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2338                    if(primalColumnSolution[iColumn]>0.5&&model.solver()->isInteger(iColumn)) {
2339                      if (comma)
2340                        fprintf(fp,",");
2341                      if (newLine)
2342                        fprintf(fp,"\n");
2343                      fprintf(fp,"%d ",iColumn);
2344                      comma=true;
2345                      newLine=false;
2346                      n++;
2347                      if (n==10) {
2348                        n=0;
2349                        newLine=true;
2350                      }
2351                    }
2352                  }
2353                  fprintf(fp,"};\n");
2354                  n=0;
2355                  comma=false;
2356                  newLine=false;
2357                  fprintf(fp,"\tdouble intSolnV[]={\n");
2358                  for ( iColumn=0;iColumn<numberColumns;iColumn++) {
2359                    if(primalColumnSolution[iColumn]>0.5&&model.solver()->isInteger(iColumn)) {
2360                      if (comma)
2361                        fprintf(fp,",");
2362                      if (newLine)
2363                        fprintf(fp,"\n");
2364                      int value = (int) (primalColumnSolution[iColumn]+0.5);
2365                      fprintf(fp,"%d. ",value);
2366                      comma=true;
2367                      newLine=false;
2368                      n++;
2369                      if (n==10) {
2370                        n=0;
2371                        newLine=true;
2372                      }
2373                    }
2374                  }
2375                  fprintf(fp,"};\n");
2376                }
2377                if (fp!=stdout)
2378                  fclose(fp);
2379              } else {
2380                std::cout<<"Unable to open file "<<fileName<<std::endl;
2381              }
2382            } else {
2383              std::cout<<"** Current model not valid"<<std::endl;
2384             
2385            }
2386            break;
2387          case SAVESOL:
2388            if (goodModel) {
2389              // get next field
2390              field = CoinReadGetString(argc,argv);
2391              if (field=="$") {
2392                field = parameters[iParam].stringValue();
2393              } else if (field=="EOL") {
2394                parameters[iParam].printString();
2395                break;
2396              } else {
2397                parameters[iParam].setStringValue(field);
2398              }
2399              std::string fileName;
2400              if (field[0]=='/'||field[0]=='\\') {
2401                fileName = field;
2402              } else if (field[0]=='~') {
2403                char * environ = getenv("HOME");
2404                if (environ) {
2405                  std::string home(environ);
2406                  field=field.erase(0,1);
2407                  fileName = home+field;
2408                } else {
2409                  fileName=field;
2410                }
2411              } else {
2412                fileName = directory+field;
2413              }
2414              saveSolution(lpSolver,fileName);
2415            } else {
2416              std::cout<<"** Current model not valid"<<std::endl;
2417             
2418            }
2419            break;
2420          default:
2421            abort();
2422          }
2423        } 
2424      } else if (!numberMatches) {
2425        std::cout<<"No match for "<<field<<" - ? for list of commands"
2426                 <<std::endl;
2427      } else if (numberMatches==1) {
2428        if (!numberQuery) {
2429          std::cout<<"Short match for "<<field<<" - completion: ";
2430          std::cout<<parameters[firstMatch].matchName()<<std::endl;
2431        } else if (numberQuery) {
2432          std::cout<<parameters[firstMatch].matchName()<<" : ";
2433          std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
2434          if (numberQuery>=2) 
2435            parameters[firstMatch].printLongHelp();
2436        }
2437      } else {
2438        if (!numberQuery) 
2439          std::cout<<"Multiple matches for "<<field<<" - possible completions:"
2440                   <<std::endl;
2441        else
2442          std::cout<<"Completions of "<<field<<":"<<std::endl;
2443        for ( iParam=0; iParam<numberParameters; iParam++ ) {
2444          int match = parameters[iParam].matches(field);
2445          if (match&&parameters[iParam].displayThis()) {
2446            std::cout<<parameters[iParam].matchName();
2447            if (numberQuery>=2) 
2448              std::cout<<" : "<<parameters[iParam].shortHelp();
2449            std::cout<<std::endl;
2450          }
2451        }
2452      }
2453    }
2454  }
2455  // By now all memory should be freed
2456#ifdef DMALLOC
2457  dmalloc_log_unfreed();
2458  dmalloc_shutdown();
2459#endif
2460  return 0;
2461}   
Note: See TracBrowser for help on using the repository browser.