source: branches/devel/Cbc/src/CbcSolver.cpp @ 590

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

new scheme doesn't like batch files

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