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

Last change on this file since 616 was 615, checked in by forrest, 12 years ago

for readLP names

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