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

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

conflict

File size: 230.9 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.setAbsoluteIncrement(4332.64);
2708                      cbcModel->addHeuristic(&heuristicFPump);
2709                      heuristicFPump.setInitialWeight(1);
2710                     
2711                      CbcRounding rounding(*cbcModel);
2712                      cbcModel->addHeuristic(&rounding);
2713                     
2714                      CbcHeuristicLocal heuristicLocal(*cbcModel);
2715                      heuristicLocal.setSearchType(1);
2716                      cbcModel->addHeuristic(&heuristicLocal);
2717                     
2718                      CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
2719                      cbcModel->addHeuristic(&heuristicGreedyCover);
2720                     
2721                      CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
2722                      cbcModel->addHeuristic(&heuristicGreedyEquality);
2723                     
2724                      CbcCompareDefault compare;
2725                      cbcModel->setNodeComparison(compare);
2726                      cbcModel->setNumberBeforeTrust(5);
2727                      cbcModel->setSpecialOptions(2);
2728                      cbcModel->messageHandler()->setLogLevel(1);
2729                      cbcModel->setMaximumCutPassesAtRoot(-100);
2730                      cbcModel->setMaximumCutPasses(1);
2731                      cbcModel->setMinimumDrop(0.05);
2732                      // For branchAndBound this may help
2733                      clpModel->defaultFactorizationFrequency();
2734                      clpModel->setDualBound(1.0001e+08);
2735                      clpModel->setPerturbation(50);
2736                      osiclpModel->setSpecialOptions(193);
2737                      osiclpModel->messageHandler()->setLogLevel(0);
2738                      osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
2739                      osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
2740                      // You can save some time by switching off message building
2741                      // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
2742                     
2743                      // Solve
2744                     
2745                      cbcModel->initialSolve();
2746                      if (clpModel->tightenPrimalBounds()!=0) {
2747                        std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
2748                        exit(1);
2749                      }
2750                      clpModel->dual();  // clean up
2751                      cbcModel->initialSolve();
2752                      cbcModel->branchAndBound();
2753                      OsiSolverLinearizedQuadratic * solver3 = dynamic_cast<OsiSolverLinearizedQuadratic *> (model2.solver());
2754                      assert (solver3);
2755                      solution = solver3->bestSolution();
2756                      double bestObjectiveValue = solver3->bestObjectiveValue();
2757                      linkSolver->setBestObjectiveValue(bestObjectiveValue);
2758                      linkSolver->setBestSolution(solution,solver3->getNumCols());
2759                      CbcHeuristicDynamic3 dynamic(model);
2760                      model.addHeuristic(&dynamic);
2761                      // if convex
2762                      if ((linkSolver->specialOptions2()&4)!=0) {
2763                        int numberColumns = coinModel->numberColumns();
2764                        // add OA cut
2765                        double offset;
2766                        double * gradient = new double [numberColumns+1];
2767                        memcpy(gradient,qp->objectiveAsObject()->gradient(qp,solution,offset,true,2),
2768                               numberColumns*sizeof(double));
2769                        double rhs = 0.0;
2770                        int * column = new int[numberColumns+1];
2771                        int n=0;
2772                        for (int i=0;i<numberColumns;i++) {
2773                          double value = gradient[i];
2774                          if (fabs(value)>1.0e-12) {
2775                            gradient[n]=value;
2776                            rhs += value*solution[i];
2777                            column[n++]=i;
2778                          }
2779                        }
2780                        gradient[n]=-1.0;
2781                        column[n++]=numberColumns;
2782                        storedAmpl.addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
2783                        delete [] gradient;
2784                        delete [] column;
2785                      }
2786                      // could do three way branching round a) continuous b) best solution
2787                      printf("obj %g\n",bestObjectiveValue);
2788                      linkSolver->initialSolve();
2789                    }
2790                  }
2791                }
2792                si->setSpecialOptions(0x40000000);
2793              }
2794#endif
2795              if (!miplib) {
2796                if (!preSolve) {
2797                  model.solver()->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry);
2798                  model.solver()->setHintParam(OsiDoPresolveInResolve,false,OsiHintTry);
2799                }
2800                model.initialSolve();
2801                OsiSolverInterface * solver = model.solver();
2802                OsiClpSolverInterface * si =
2803                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
2804                ClpSimplex * clpSolver = si->getModelPtr();
2805                clpSolver->setSpecialOptions(clpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
2806                if (!complicatedInteger&&clpSolver->tightenPrimalBounds()!=0) {
2807                  std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
2808                  exit(1);
2809                }
2810                if (clpSolver->dualBound()==1.0e10) {
2811                  // user did not set - so modify
2812                  // get largest scaled away from bound
2813                  double largest=1.0e-12;
2814                  int numberRows = clpSolver->numberRows();
2815                  const double * rowPrimal = clpSolver->primalRowSolution();
2816                  const double * rowLower = clpSolver->rowLower();
2817                  const double * rowUpper = clpSolver->rowUpper();
2818                  const double * rowScale = clpSolver->rowScale();
2819                  int iRow;
2820                  for (iRow=0;iRow<numberRows;iRow++) {
2821                    double value = rowPrimal[iRow];
2822                    double above = value-rowLower[iRow];
2823                    double below = rowUpper[iRow]-value;
2824                    if (rowScale) {
2825                      double multiplier = rowScale[iRow];
2826                      above *= multiplier;
2827                      below *= multiplier;
2828                    }
2829                    if (above<1.0e12)
2830                      largest = CoinMax(largest,above);
2831                    if (below<1.0e12)
2832                      largest = CoinMax(largest,below);
2833                  }
2834                 
2835                  int numberColumns = clpSolver->numberColumns();
2836                  const double * columnPrimal = clpSolver->primalColumnSolution();
2837                  const double * columnLower = clpSolver->columnLower();
2838                  const double * columnUpper = clpSolver->columnUpper();
2839                  const double * columnScale = clpSolver->columnScale();
2840                  int iColumn;
2841                  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2842                    double value = columnPrimal[iColumn];
2843                    double above = value-columnLower[iColumn];
2844                    double below = columnUpper[iColumn]-value;
2845                    if (columnScale) {
2846                      double multiplier = 1.0/columnScale[iColumn];
2847                      above *= multiplier;
2848                      below *= multiplier;
2849                    }
2850                    if (above<1.0e12)
2851                      largest = CoinMax(largest,above);
2852                    if (below<1.0e12)
2853                      largest = CoinMax(largest,below);
2854                  }
2855                  //if (!noPrinting)
2856                  //std::cout<<"Largest (scaled) away from bound "<<largest<<std::endl;
2857                  clpSolver->setDualBound(CoinMax(1.0001e8,CoinMin(1000.0*largest,1.00001e10)));
2858                }
2859                clpSolver->dual();  // clean up
2860              }
2861              // If user made settings then use them
2862              if (!defaultSettings) {
2863                OsiSolverInterface * solver = model.solver();
2864                if (!doScaling)
2865                  solver->setHintParam(OsiDoScale,false,OsiHintTry);
2866                OsiClpSolverInterface * si =
2867                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
2868                assert (si != NULL);
2869                // get clp itself
2870                ClpSimplex * modelC = si->getModelPtr();
2871                //if (modelC->tightenPrimalBounds()!=0) {
2872                //std::cout<<"Problem is infeasible!"<<std::endl;
2873                //break;
2874                //}
2875                // bounds based on continuous
2876                if (tightenFactor&&!complicatedInteger) {
2877                  if (modelC->tightenPrimalBounds(tightenFactor)!=0) {
2878                    std::cout<<"Problem is infeasible!"<<std::endl;
2879                    break;
2880                  }
2881                }
2882                if (djFix<1.0e20) {
2883                  // do some fixing
2884                  int numberColumns = modelC->numberColumns();
2885                  int i;
2886                  const char * type = modelC->integerInformation();
2887                  double * lower = modelC->columnLower();
2888                  double * upper = modelC->columnUpper();
2889                  double * solution = modelC->primalColumnSolution();
2890                  double * dj = modelC->dualColumnSolution();
2891                  int numberFixed=0;
2892                  for (i=0;i<numberColumns;i++) {
2893                    if (type[i]) {
2894                      double value = solution[i];
2895                      if (value<lower[i]+1.0e-5&&dj[i]>djFix) {
2896                        solution[i]=lower[i];
2897                        upper[i]=lower[i];
2898                        numberFixed++;
2899                      } else if (value>upper[i]-1.0e-5&&dj[i]<-djFix) {
2900                        solution[i]=upper[i];
2901                        lower[i]=upper[i];
2902                        numberFixed++;
2903                      }
2904                    }
2905                  }
2906                  printf("%d columns fixed\n",numberFixed);
2907                }
2908              }
2909              // See if we want preprocessing
2910              OsiSolverInterface * saveSolver=NULL;
2911              CglPreProcess process;
2912              delete babModel;
2913              babModel = new CbcModel(model);
2914              OsiSolverInterface * solver3 = clpSolver->clone();
2915              babModel->assignSolver(solver3);
2916              OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
2917              int numberChanged=0;
2918              if (clpSolver2->messageHandler()->logLevel())
2919                clpSolver2->messageHandler()->setLogLevel(1);
2920              if (logLevel>-1)
2921                clpSolver2->messageHandler()->setLogLevel(logLevel);
2922              lpSolver = clpSolver2->getModelPtr();
2923              if (lpSolver->factorizationFrequency()==200&&!miplib) {
2924                // User did not touch preset
2925                int numberRows = lpSolver->numberRows();
2926                const int cutoff1=10000;
2927                const int cutoff2=100000;
2928                const int base=75;
2929                const int freq0 = 50;
2930                const int freq1=200;
2931                const int freq2=400;
2932                const int maximum=1000;
2933                int frequency;
2934                if (numberRows<cutoff1)
2935                  frequency=base+numberRows/freq0;
2936                else if (numberRows<cutoff2)
2937                  frequency=base+cutoff1/freq0 + (numberRows-cutoff1)/freq1;
2938                else
2939                  frequency=base+cutoff1/freq0 + (cutoff2-cutoff1)/freq1 + (numberRows-cutoff2)/freq2;
2940                lpSolver->setFactorizationFrequency(CoinMin(maximum,frequency));
2941              }
2942              time2 = CoinCpuTime();
2943              totalTime += time2-time1;
2944              time1 = time2;
2945              double timeLeft = babModel->getMaximumSeconds();
2946              int numberOriginalColumns = babModel->solver()->getNumCols();
2947              if (preProcess==7) {
2948                // use strategy instead
2949                preProcess=0;
2950                useStrategy=true;
2951#ifdef COIN_HAS_LINK
2952                // empty out any cuts
2953                if (storedAmpl.sizeRowCuts()) {
2954                  printf("Emptying ampl stored cuts as internal preprocessing\n");
2955                  CglStored temp;
2956                  storedAmpl=temp;
2957                }
2958#endif
2959              }
2960              if (preProcess&&type==BAB) {
2961                // See if sos from mps file
2962                if (numberSOS==0&&clpSolver->numberSOS()&&doSOS) {
2963                  // SOS
2964                  numberSOS = clpSolver->numberSOS();
2965                  const CoinSet * setInfo = clpSolver->setInfo();
2966                  sosStart = new int [numberSOS+1];
2967                  sosType = new char [numberSOS];
2968                  int i;
2969                  int nTotal=0;
2970                  sosStart[0]=0;
2971                  for ( i=0;i<numberSOS;i++) {
2972                    int type = setInfo[i].setType();
2973                    int n=setInfo[i].numberEntries();
2974                    sosType[i]=type;
2975                    nTotal += n;
2976                    sosStart[i+1] = nTotal;
2977                  }
2978                  sosIndices = new int[nTotal];
2979                  sosReference = new double [nTotal];
2980                  for (i=0;i<numberSOS;i++) {
2981                    int n=setInfo[i].numberEntries();
2982                    const int * which = setInfo[i].which();
2983                    const double * weights = setInfo[i].weights();
2984                    int base = sosStart[i];
2985                    for (int j=0;j<n;j++) {
2986                      int k=which[j];
2987                      sosIndices[j+base]=k;
2988                      sosReference[j+base] = weights ? weights[j] : (double) j;
2989                    }
2990                  }
2991                }
2992                saveSolver=babModel->solver()->clone();
2993                /* Do not try and produce equality cliques and
2994                   do up to 10 passes */
2995                OsiSolverInterface * solver2;
2996                {
2997                  // Tell solver we are in Branch and Cut
2998                  saveSolver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo) ;
2999                  // Default set of cut generators
3000                  CglProbing generator1;
3001                  generator1.setUsingObjective(1);
3002                  generator1.setMaxPass(3);
3003                  generator1.setMaxProbeRoot(saveSolver->getNumCols());
3004                  generator1.setMaxElements(100);
3005                  generator1.setMaxLookRoot(50);
3006                  generator1.setRowCuts(3);
3007                  // Add in generators
3008                  process.addCutGenerator(&generator1);
3009                  int translate[]={9999,0,0,-1,2,3,-2};
3010                  process.messageHandler()->setLogLevel(babModel->logLevel());
3011#ifdef COIN_HAS_ASL
3012                  if (info.numberSos&&doSOS&&usingAmpl) {
3013                    // SOS
3014                    numberSOS = info.numberSos;
3015                    sosStart = info.sosStart;
3016                    sosIndices = info.sosIndices;
3017                  }
3018#endif
3019                  if (numberSOS&&doSOS) {
3020                    // SOS
3021                    int numberColumns = saveSolver->getNumCols();
3022                    char * prohibited = new char[numberColumns];
3023                    memset(prohibited,0,numberColumns);
3024                    int n=sosStart[numberSOS];
3025                    for (int i=0;i<n;i++) {
3026                      int iColumn = sosIndices[i];
3027                      prohibited[iColumn]=1;
3028                    }
3029                    process.passInProhibited(prohibited,numberColumns);
3030                    delete [] prohibited;
3031                  }
3032                  int numberPasses = 10;
3033                  if (tunePreProcess>=1000000) {
3034                    numberPasses = (tunePreProcess/1000000)-1;
3035                    tunePreProcess = tunePreProcess % 1000000;
3036                  } else if (tunePreProcess>=1000) {
3037                    numberPasses = (tunePreProcess/1000)-1;
3038                    tunePreProcess = tunePreProcess % 1000;
3039                  }
3040                  if (doSprint>0) {
3041                    // Sprint for primal solves
3042                    ClpSolve::SolveType method = ClpSolve::usePrimalorSprint;
3043                    ClpSolve::PresolveType presolveType = ClpSolve::presolveOff;
3044                    int numberPasses = 5;
3045                    int options[] = {0,3,0,0,0,0};
3046                    int extraInfo[] = {-1,20,-1,-1,-1,-1};
3047                    extraInfo[1]=doSprint;
3048                    int independentOptions[] = {0,0,3};
3049                    ClpSolve clpSolve(method,presolveType,numberPasses,
3050                                      options,extraInfo,independentOptions);
3051                    // say use in OsiClp
3052                    clpSolve.setSpecialOption(6,1);
3053                    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (saveSolver);
3054                    osiclp->setSolveOptions(clpSolve);
3055                    osiclp->setHintParam(OsiDoDualInResolve,false);
3056                    // switch off row copy
3057                    osiclp->getModelPtr()->setSpecialOptions(osiclp->getModelPtr()->specialOptions()|256);
3058                    osiclp->getModelPtr()->setInfeasibilityCost(1.0e11);
3059                  }
3060                  solver2 = process.preProcessNonDefault(*saveSolver,translate[preProcess],numberPasses,
3061                                                         tunePreProcess);
3062                  // Tell solver we are not in Branch and Cut
3063                  saveSolver->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
3064                  if (solver2)
3065                    solver2->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo) ;
3066                }
3067#ifdef COIN_HAS_ASL
3068                if (!solver2&&usingAmpl) {
3069                  // infeasible
3070                  info.problemStatus=1;
3071                  info.objValue = 1.0e100;
3072                  sprintf(info.buffer,"infeasible/unbounded by pre-processing");
3073                  info.primalSolution=NULL;
3074                  info.dualSolution=NULL;
3075                  break;
3076                }
3077#endif
3078                if (!noPrinting) {
3079                  if (!solver2) {
3080                    printf("Pre-processing says infeasible or unbounded\n");
3081                    break;
3082                  } else {
3083                    printf("processed model has %d rows, %d columns and %d elements\n",
3084                           solver2->getNumRows(),solver2->getNumCols(),solver2->getNumElements());
3085                  }
3086                }
3087                //solver2->resolve();
3088                if (preProcess==2) {
3089                  OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (solver2);
3090                  ClpSimplex * lpSolver = clpSolver2->getModelPtr();
3091                  lpSolver->writeMps("presolved.mps",0,1,lpSolver->optimizationDirection());
3092                  printf("Preprocessed model (minimization) on presolved.mps\n");
3093                }
3094                // we have to keep solver2 so pass clone
3095                solver2 = solver2->clone();
3096                babModel->assignSolver(solver2);
3097                babModel->initialSolve();
3098                babModel->setMaximumSeconds(timeLeft-(CoinCpuTime()-time1));
3099              }
3100              // now tighten bounds
3101              if (!miplib) {
3102                OsiClpSolverInterface * si =
3103                  dynamic_cast<OsiClpSolverInterface *>(babModel->solver()) ;
3104                assert (si != NULL);
3105                // get clp itself
3106                ClpSimplex * modelC = si->getModelPtr();
3107                if (noPrinting)
3108                  modelC->setLogLevel(0);
3109                if (!complicatedInteger&&modelC->tightenPrimalBounds()!=0) {
3110                  std::cout<<"Problem is infeasible!"<<std::endl;
3111                  break;
3112                }
3113                modelC->dual();
3114              }
3115#if 0
3116              numberDebugValues=599;
3117              debugValues = new double[numberDebugValues];
3118              CoinZeroN(debugValues,numberDebugValues);
3119              debugValues[3]=1.0;
3120              debugValues[6]=25.0;
3121              debugValues[9]=4.0;
3122              debugValues[26]=4.0;
3123              debugValues[27]=6.0;
3124              debugValues[35]=8.0;
3125              debugValues[53]=21.0;
3126              debugValues[56]=4.0;
3127#endif
3128              if (debugValues) {
3129                // for debug
3130                std::string problemName ;
3131                babModel->solver()->getStrParam(OsiProbName,problemName) ;
3132                babModel->solver()->activateRowCutDebugger(problemName.c_str()) ;
3133                twomirGen.probname_=strdup(problemName.c_str());
3134                // checking seems odd
3135                //redsplitGen.set_given_optsol(babModel->solver()->getRowCutDebuggerAlways()->optimalSolution(),
3136                //                         babModel->getNumCols());
3137              }
3138              int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
3139#ifdef COIN_HAS_LINK
3140              // If linked then see if expansion wanted
3141              {
3142                OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
3143                if (solver3) {
3144                  int options = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()/10000;
3145                  if (options) {
3146                    /*
3147                      1 - force mini branch and bound
3148                      2 - set priorities high on continuous
3149                      4 - try adding OA cuts
3150                      8 - try doing quadratic linearization
3151                      16 - try expanding knapsacks
3152                    */
3153                    if ((options&16)) {
3154                      int numberColumns = saveCoinModel.numberColumns();
3155                      int numberRows = saveCoinModel.numberRows();
3156                      whichColumn = new int[numberColumns];
3157                      knapsackStart=new int[numberRows+1];
3158                      knapsackRow=new int[numberRows];
3159                      numberKnapsack=10000;
3160                      int extra1 = parameters[whichParam(EXTRA1,numberParameters,parameters)].intValue();
3161                      int extra2 = parameters[whichParam(EXTRA2,numberParameters,parameters)].intValue();
3162                      int logLevel = parameters[log].intValue();
3163                      OsiSolverInterface * solver = expandKnapsack(saveCoinModel,whichColumn,knapsackStart,
3164                                                                   knapsackRow,numberKnapsack,
3165                                                                   storedAmpl,logLevel,extra1,extra2,
3166                                                                   saveTightenedModel);
3167                      if (solver) {
3168                        clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3169                        assert (clpSolver);
3170                        lpSolver = clpSolver->getModelPtr();
3171                        babModel->assignSolver(solver);
3172                        testOsiOptions=0;
3173                        // Priorities already done
3174                        free(info.priorities);
3175                        info.priorities=NULL;
3176                      } else {
3177                        numberKnapsack=0;
3178                        delete [] whichColumn;
3179                        delete [] knapsackStart;
3180                        delete [] knapsackRow;
3181                        whichColumn=NULL;
3182                        knapsackStart=NULL;
3183                        knapsackRow=NULL;
3184                      }
3185                    }
3186                  }
3187                }
3188              }
3189#endif
3190              if (useCosts&&testOsiOptions<0) {
3191                int numberColumns = babModel->getNumCols();
3192                int * sort = new int[numberColumns];
3193                double * dsort = new double[numberColumns];
3194                int * priority = new int [numberColumns];
3195                const double * objective = babModel->getObjCoefficients();
3196                const double * lower = babModel->getColLower() ;
3197                const double * upper = babModel->getColUpper() ;
3198                const CoinPackedMatrix * matrix = babModel->solver()->getMatrixByCol();
3199                const int * columnLength = matrix->getVectorLengths();
3200                int iColumn;
3201                int n=0;
3202                for (iColumn=0;iColumn<numberColumns;iColumn++) {
3203                  if (babModel->isInteger(iColumn)) {
3204                    sort[n]=n;
3205                    if (useCosts==1)
3206                      dsort[n++]=-fabs(objective[iColumn]);
3207                    else if (useCosts==2)
3208                      dsort[n++]=iColumn;
3209                    else if (useCosts==3)
3210                      dsort[n++]=upper[iColumn]-lower[iColumn];
3211                    else if (useCosts==4)
3212                      dsort[n++]=-(upper[iColumn]-lower[iColumn]);
3213                    else if (useCosts==5)
3214                      dsort[n++]=-columnLength[iColumn];
3215                  }
3216                }
3217                CoinSort_2(dsort,dsort+n,sort);
3218                int level=0;
3219                double last = -1.0e100;
3220                for (int i=0;i<n;i++) {
3221                  int iPut=sort[i];
3222                  if (dsort[i]!=last) {
3223                    level++;
3224                    last=dsort[i];
3225                  }
3226                  priority[iPut]=level;
3227                }
3228                babModel->passInPriorities( priority,false);
3229                delete [] priority;
3230                delete [] sort;
3231                delete [] dsort;
3232              }
3233              // FPump done first as it only works if no solution
3234              CbcHeuristicFPump heuristic4(*babModel);
3235              if (useFpump) {
3236                heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
3237                int pumpTune=parameters[whichParam(FPUMPTUNE,numberParameters,parameters)].intValue();
3238                if (pumpTune>0) {
3239                  /*
3240                    >=10000000 for using obj
3241                    >=1000000 use as accumulate switch
3242                    >=1000 use index+1 as number of large loops
3243                    >=100 use 0.05 objvalue as increment
3244                    >=10 use +0.1 objvalue for cutoff (add)
3245                    1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
3246                    4 and static continuous, 5 as 3 but no internal integers
3247                    6 as 3 but all slack basis!
3248                  */
3249                  int logLevel = parameters[log].intValue();
3250                  double value = babModel->solver()->getObjSense()*babModel->solver()->getObjValue();
3251                  int w = pumpTune/10;
3252                  int c = w % 10;
3253                  w /= 10;
3254                  int i = w % 10;
3255                  w /= 10;
3256                  int r = w;
3257                  int accumulate = r/1000;
3258                  r -= 1000*accumulate;
3259                  if (accumulate>=10) {
3260                    int which = accumulate/10;
3261                    accumulate -= 10*which;
3262                    which--;
3263                    // weights and factors
3264                    double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
3265                    double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
3266                    heuristic4.setInitialWeight(weight[which]);
3267                    heuristic4.setWeightFactor(factor[which]);
3268                  }
3269                  // fake cutoff
3270                  if (logLevel>1)
3271                    printf("Setting ");
3272                  if (c) {
3273                    double cutoff;
3274                    babModel->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
3275                    cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
3276                    heuristic4.setFakeCutoff(cutoff);
3277                    if (logLevel>1)
3278                      printf("fake cutoff of %g ",cutoff);
3279                  }
3280                  if (i||r) {
3281                    // also set increment
3282                    heuristic4.setAbsoluteIncrement((0.01*i+0.005)*(fabs(value)+1.0e-12));
3283                    heuristic4.setAccumulate(accumulate);
3284                    heuristic4.setMaximumRetries(r+1);
3285                    if (logLevel>1) {
3286                      if (i) 
3287                        printf("increment of %g ",heuristic4.absoluteIncrement());
3288                      if (accumulate) 
3289                        printf("accumulate of %d ",accumulate);
3290                      printf("%d retries ",r+2);
3291                    }
3292                  }
3293                  pumpTune = pumpTune%100;
3294                  if (logLevel>1)
3295                    printf("and setting when to %d\n",pumpTune+10);
3296                  if (pumpTune==6)
3297                    pumpTune =13;
3298                  heuristic4.setWhen(pumpTune+10);
3299                }
3300                babModel->addHeuristic(&heuristic4);
3301              }
3302              if (!miplib) {
3303                CbcRounding heuristic1(*babModel);
3304                if (useRounding)
3305                  babModel->addHeuristic(&heuristic1) ;
3306                CbcHeuristicLocal heuristic2(*babModel);
3307                heuristic2.setSearchType(1);
3308                if (useCombine)
3309                  babModel->addHeuristic(&heuristic2);
3310                CbcHeuristicGreedyCover heuristic3(*babModel);
3311                CbcHeuristicGreedyEquality heuristic3a(*babModel);
3312                if (useGreedy) {
3313                  babModel->addHeuristic(&heuristic3);
3314                  babModel->addHeuristic(&heuristic3a);
3315                }
3316                if (useLocalTree) {
3317                  CbcTreeLocal localTree(babModel,NULL,10,0,0,10000,2000);
3318                  babModel->passInTreeHandler(localTree);
3319                }
3320              }
3321              CbcHeuristicRINS heuristic5(*babModel);
3322              if (useRINS)
3323                babModel->addHeuristic(&heuristic5) ;
3324              if (type==MIPLIB) {
3325                if (babModel->numberStrong()==5&&babModel->numberBeforeTrust()==5) 
3326                  babModel->setNumberBeforeTrust(50);
3327              }
3328              // add cut generators if wanted
3329              int switches[20];
3330              int numberGenerators=0;
3331              int translate[]={-100,-1,-99,-98,1,1,1,1};
3332              if (probingAction) {
3333                if (probingAction==5||probingAction==7)
3334                  probingGen.setRowCuts(-3); // strengthening etc just at root
3335                if (probingAction==6||probingAction==7) {
3336                  // Number of unsatisfied variables to look at
3337                  probingGen.setMaxProbe(1000);
3338                  probingGen.setMaxProbeRoot(1000);
3339                  // How far to follow the consequences
3340                  probingGen.setMaxLook(50);
3341                  probingGen.setMaxLookRoot(50);
3342                }
3343                babModel->addCutGenerator(&probingGen,translate[probingAction],"Probing");
3344                switches[numberGenerators++]=0;
3345              }
3346              if (gomoryAction&&(complicatedInteger!=1||gomoryAction==1)) {
3347                babModel->addCutGenerator(&gomoryGen,translate[gomoryAction],"Gomory");
3348                switches[numberGenerators++]=-1;
3349              }
3350              if (knapsackAction) {
3351                babModel->addCutGenerator(&knapsackGen,translate[knapsackAction],"Knapsack");
3352                switches[numberGenerators++]=0;
3353              }
3354              if (redsplitAction&&!complicatedInteger) {
3355                babModel->addCutGenerator(&redsplitGen,translate[redsplitAction],"Reduce-and-split");
3356                switches[numberGenerators++]=1;
3357              }
3358              if (cliqueAction) {
3359                babModel->addCutGenerator(&cliqueGen,translate[cliqueAction],"Clique");
3360                switches[numberGenerators++]=0;
3361              }
3362              if (mixedAction) {
3363                babModel->addCutGenerator(&mixedGen,translate[mixedAction],"MixedIntegerRounding2");
3364                switches[numberGenerators++]=-1;
3365              }
3366              if (flowAction) {
3367                babModel->addCutGenerator(&flowGen,translate[flowAction],"FlowCover");
3368                switches[numberGenerators++]=1;
3369              }
3370              if (twomirAction&&!complicatedInteger) {
3371                babModel->addCutGenerator(&twomirGen,translate[twomirAction],"TwoMirCuts");
3372                switches[numberGenerators++]=1;
3373              }
3374              if (landpAction) {
3375                babModel->addCutGenerator(&landpGen,translate[landpAction],"LiftAndProject");
3376                switches[numberGenerators++]=1;
3377              }
3378              if (storedCuts) 
3379                babModel->setSpecialOptions(babModel->specialOptions()|64);
3380              // Say we want timings
3381              numberGenerators = babModel->numberCutGenerators();
3382              int iGenerator;
3383              int cutDepth=
3384                parameters[whichParam(CUTDEPTH,numberParameters,parameters)].intValue();
3385              for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
3386                CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
3387                int howOften = generator->howOften();
3388                if (howOften==-98||howOften==-99) 
3389                  generator->setSwitchOffIfLessThan(switches[iGenerator]);
3390                generator->setTiming(true);
3391                if (cutDepth>=0)
3392                  generator->setWhatDepth(cutDepth) ;
3393              }
3394              // Could tune more
3395              if (!miplib) {
3396                babModel->setMinimumDrop(min(5.0e-2,
3397                                             fabs(babModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
3398                if (cutPass==-1234567) {
3399                  if (babModel->getNumCols()<500)
3400                    babModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
3401                  else if (babModel->getNumCols()<5000)
3402                    babModel->setMaximumCutPassesAtRoot(100); // use minimum drop
3403                  else
3404                    babModel->setMaximumCutPassesAtRoot(20);
3405                } else {
3406                  babModel->setMaximumCutPassesAtRoot(cutPass);
3407                }
3408                babModel->setMaximumCutPasses(1);
3409              }
3410              // Do more strong branching if small
3411              //if (babModel->getNumCols()<5000)
3412              //babModel->setNumberStrong(20);
3413              // Switch off strong branching if wanted
3414              //if (babModel->getNumCols()>10*babModel->getNumRows())
3415              //babModel->setNumberStrong(0);
3416              if (!noPrinting) {
3417                int iLevel = parameters[log].intValue();
3418                if (iLevel<0) {
3419                  babModel->setPrintingMode(1);
3420                  iLevel = -iLevel;
3421                }
3422                babModel->messageHandler()->setLogLevel(iLevel);
3423                if (babModel->getNumCols()>2000||babModel->getNumRows()>1500||
3424                    babModel->messageHandler()->logLevel()>1)
3425                  babModel->setPrintFrequency(100);
3426              }
3427             
3428              babModel->solver()->setIntParam(OsiMaxNumIterationHotStart,
3429                    parameters[whichParam(MAXHOTITS,numberParameters,parameters)].intValue());
3430              OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
3431              // go faster stripes
3432              if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
3433                osiclp->setupForRepeatedUse(2,parameters[slog].intValue());
3434              } else {
3435                osiclp->setupForRepeatedUse(0,parameters[slog].intValue());
3436              }
3437              double increment=babModel->getCutoffIncrement();;
3438              int * changed = NULL;
3439              if (!miplib&&increment==normalIncrement)
3440                changed=analyze( osiclp,numberChanged,increment,false);
3441              if (debugValues) {
3442                if (numberDebugValues==babModel->getNumCols()) {
3443                  // for debug
3444                  babModel->solver()->activateRowCutDebugger(debugValues) ;
3445                } else {
3446                  printf("debug file has incorrect number of columns\n");
3447                }
3448              }
3449              babModel->setCutoffIncrement(CoinMax(babModel->getCutoffIncrement(),increment));
3450              // Turn this off if you get problems
3451              // Used to be automatically set
3452              int mipOptions = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()%10000;
3453              if (mipOptions!=(1))
3454                printf("mip options %d\n",mipOptions);
3455              osiclp->setSpecialOptions(mipOptions);
3456              if (gapRatio < 1.0e100) {
3457                double value = babModel->solver()->getObjValue() ;
3458                double value2 = gapRatio*(1.0e-5+fabs(value)) ;
3459                babModel->setAllowableGap(value2) ;
3460                std::cout << "Continuous " << value
3461                          << ", so allowable gap set to "
3462                          << value2 << std::endl ;
3463              }
3464              // probably faster to use a basis to get integer solutions
3465              babModel->setSpecialOptions(babModel->specialOptions()|2);
3466              currentBranchModel = babModel;
3467              OsiSolverInterface * strengthenedModel=NULL;
3468              if (type==BAB||type==MIPLIB) {
3469                int moreMipOptions = parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].intValue();
3470                if (moreMipOptions>=0) {
3471                  printf("more mip options %d\n",moreMipOptions);
3472                  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
3473                  if (moreMipOptions==10000) {
3474                    // test memory saving
3475                    moreMipOptions -= 10000;
3476                    ClpSimplex * lpSolver = osiclp->getModelPtr();
3477                    lpSolver->setPersistenceFlag(1);
3478                    // switch off row copy if few rows
3479                    if (lpSolver->numberRows()<150)
3480                      lpSolver->setSpecialOptions(lpSolver->specialOptions()|256);
3481                  }
3482                  if (((moreMipOptions+1)%1000000)!=0)
3483                    babModel->setSearchStrategy(moreMipOptions%1000000);
3484                  // go faster stripes
3485                  if( moreMipOptions >=999999) {
3486                    if (osiclp) {
3487                      int save = osiclp->specialOptions();
3488                      osiclp->setupForRepeatedUse(2,0);
3489                      osiclp->setSpecialOptions(save|osiclp->specialOptions());
3490                    }
3491                  } 
3492                }
3493              }
3494              if (type==BAB) {
3495#ifdef COIN_HAS_ASL
3496                if (usingAmpl) {
3497                  priorities=info.priorities;
3498                  branchDirection=info.branchDirection;
3499                  pseudoDown=info.pseudoDown;
3500                  pseudoUp=info.pseudoUp;
3501                  solutionIn=info.primalSolution;
3502                  prioritiesIn = info.priorities;
3503                  if (info.numberSos&&doSOS) {
3504                    // SOS
3505                    numberSOS = info.numberSos;
3506                    sosStart = info.sosStart;
3507                    sosIndices = info.sosIndices;
3508                    sosType = info.sosType;
3509                    sosReference = info.sosReference;
3510                    sosPriority = info.sosPriority;
3511                  }
3512                }
3513#endif               
3514                const int * originalColumns = preProcess ? process.originalColumns() : NULL;
3515                if (solutionIn&&useSolution) {
3516                  if (preProcess) {
3517                    int numberColumns = babModel->getNumCols();
3518                    // extend arrays in case SOS
3519                    int n = originalColumns[numberColumns-1]+1;
3520                    int nSmaller = CoinMin(n,numberOriginalColumns);
3521                    double * solutionIn2 = new double [n];
3522                    int * prioritiesIn2 = new int[n];
3523                    int i;
3524                    for (i=0;i<nSmaller;i++) {
3525                      solutionIn2[i]=solutionIn[i];
3526                      prioritiesIn2[i]=prioritiesIn[i];
3527                    }
3528                    for (;i<n;i++) {
3529                      solutionIn2[i]=0.0;
3530                      prioritiesIn2[i]=1000000;
3531                    }
3532                    int iLast=-1;
3533                    for (i=0;i<numberColumns;i++) {
3534                      int iColumn = originalColumns[i];
3535                      assert (iColumn>iLast);
3536                      iLast=iColumn;
3537                      solutionIn2[i]=solutionIn2[iColumn];
3538                      if (prioritiesIn)
3539                        prioritiesIn2[i]=prioritiesIn2[iColumn];
3540                    }
3541                    babModel->setHotstartSolution(solutionIn2,prioritiesIn2);
3542                    delete [] solutionIn2;
3543                    delete [] prioritiesIn2;
3544                  } else {
3545                    babModel->setHotstartSolution(solutionIn,prioritiesIn);
3546                  }
3547                }
3548                OsiSolverInterface * testOsiSolver= (testOsiOptions>=0) ? babModel->solver() : NULL;
3549                if (!testOsiSolver) {
3550                  // *************************************************************
3551                  // CbcObjects
3552                  if (preProcess&&process.numberSOS()) {
3553                    int numberSOS = process.numberSOS();
3554                    int numberIntegers = babModel->numberIntegers();
3555                    /* model may not have created objects
3556                       If none then create
3557                    */
3558                    if (!numberIntegers||!babModel->numberObjects()) {
3559                      int type = (pseudoUp) ? 1 : 0;
3560                      babModel->findIntegers(true,type);
3561                      numberIntegers = babModel->numberIntegers();
3562                    }
3563                    OsiObject ** oldObjects = babModel->objects();
3564                    // Do sets and priorities
3565                    OsiObject ** objects = new OsiObject * [numberSOS];
3566                    // set old objects to have low priority
3567                    int numberOldObjects = babModel->numberObjects();
3568                    int numberColumns = babModel->getNumCols();
3569                    for (int iObj = 0;iObj<numberOldObjects;iObj++) {
3570                      oldObjects[iObj]->setPriority(numberColumns+1);
3571                      int iColumn = oldObjects[iObj]->columnNumber();
3572                      assert (iColumn>=0);
3573                      if (iColumn>=numberOriginalColumns)
3574                        continue;
3575                      if (originalColumns)
3576                        iColumn = originalColumns[iColumn];
3577                      if (branchDirection) {
3578                        CbcSimpleInteger * obj =
3579                          dynamic_cast <CbcSimpleInteger *>(oldObjects[iObj]) ;
3580                        if (obj) { 
3581                          obj->setPreferredWay(branchDirection[iColumn]);
3582                        } else {
3583                          CbcObject * obj =
3584                            dynamic_cast <CbcObject *>(oldObjects[iObj]) ;
3585                          assert (obj);
3586                          obj->setPreferredWay(branchDirection[iColumn]);
3587                        }
3588                      }
3589                      if (pseudoUp) {
3590                        CbcSimpleIntegerPseudoCost * obj1a =
3591                          dynamic_cast <CbcSimpleIntegerPseudoCost *>(oldObjects[iObj]) ;
3592                        assert (obj1a);
3593                        if (pseudoDown[iColumn]>0.0)
3594                          obj1a->setDownPseudoCost(pseudoDown[iColumn]);
3595                        if (pseudoUp[iColumn]>0.0)
3596                          obj1a->setUpPseudoCost(pseudoUp[iColumn]);
3597                      }
3598                    }
3599                    const int * starts = process.startSOS();
3600                    const int * which = process.whichSOS();
3601                    const int * type = process.typeSOS();
3602                    const double * weight = process.weightSOS();
3603                    int iSOS;
3604                    for (iSOS =0;iSOS<numberSOS;iSOS++) {
3605                      int iStart = starts[iSOS];
3606                      int n=starts[iSOS+1]-iStart;
3607                      objects[iSOS] = new CbcSOS(babModel,n,which+iStart,weight+iStart,
3608                                                 iSOS,type[iSOS]);
3609                      // branch on long sets first
3610                      objects[iSOS]->setPriority(numberColumns-n);
3611                    }
3612                    babModel->addObjects(numberSOS,objects);
3613                    for (iSOS=0;iSOS<numberSOS;iSOS++)
3614                      delete objects[iSOS];
3615                    delete [] objects;
3616                  } else if (priorities||branchDirection||pseudoDown||pseudoUp||numberSOS) {
3617                    // do anyway for priorities etc
3618                    int numberIntegers = babModel->numberIntegers();
3619                    /* model may not have created objects
3620                       If none then create
3621                    */
3622                    if (!numberIntegers||!babModel->numberObjects()) {
3623                      int type = (pseudoUp) ? 1 : 0;
3624                      babModel->findIntegers(true,type);
3625                    }
3626                    if (numberSOS) {
3627                      // Do sets and priorities
3628                      OsiObject ** objects = new OsiObject * [numberSOS];
3629                      int iSOS;
3630                      if (originalColumns) {
3631                        // redo sequence numbers
3632                        int numberColumns = babModel->getNumCols();
3633                        int nOld = originalColumns[numberColumns-1]+1;
3634                        int * back = new int[nOld];
3635                        int i;
3636                        for (i=0;i<nOld;i++)
3637                          back[i]=-1;
3638                        for (i=0;i<numberColumns;i++)
3639                          back[originalColumns[i]]=i;
3640                        // Really need better checks
3641                        int nMissing=0;
3642                        int n=sosStart[numberSOS];
3643                        for (i=0;i<n;i++) {
3644                          int iColumn = sosIndices[i];
3645                          int jColumn = back[iColumn];
3646                          if (jColumn>=0) 
3647                            sosIndices[i] = jColumn;
3648                          else 
3649                            nMissing++;
3650                        }
3651                        delete [] back;
3652                        if (nMissing)
3653                          printf("%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
3654                      }
3655                      for (iSOS =0;iSOS<numberSOS;iSOS++) {
3656                        int iStart = sosStart[iSOS];
3657                        int n=sosStart[iSOS+1]-iStart;
3658                        objects[iSOS] = new CbcSOS(babModel,n,sosIndices+iStart,sosReference+iStart,
3659                                                   iSOS,sosType[iSOS]);
3660                        if (sosPriority)
3661                          objects[iSOS]->setPriority(sosPriority[iSOS]);
3662                        else if (!prioritiesIn)
3663                          objects[iSOS]->setPriority(10);  // rather than 1000
3664                      }
3665                      // delete any existing SOS objects
3666                      int numberObjects=babModel->numberObjects();
3667                      OsiObject ** oldObjects=babModel->objects();
3668                      int nNew=0;
3669                      for (int i=0;i<numberObjects;i++) {
3670                        OsiObject * objThis = oldObjects[i];
3671                        CbcSOS * obj1 =
3672                          dynamic_cast <CbcSOS *>(objThis) ;
3673                        OsiSOS * obj2 =
3674                          dynamic_cast <OsiSOS *>(objThis) ;
3675                        if (!obj1&&!obj2) {
3676                          oldObjects[nNew++]=objThis;
3677                        } else {
3678                          delete objThis;
3679                        }
3680                      }
3681                      babModel->setNumberObjects(nNew);
3682                      babModel->addObjects(numberSOS,objects);
3683                      for (iSOS=0;iSOS<numberSOS;iSOS++)
3684                        delete objects[iSOS];
3685                      delete [] objects;
3686                    }
3687                  }
3688                  OsiObject ** objects = babModel->objects();
3689                  int numberObjects = babModel->numberObjects();
3690                  for (int iObj = 0;iObj<numberObjects;iObj++) {
3691                    // skip sos
3692                    CbcSOS * objSOS =
3693                      dynamic_cast <CbcSOS *>(objects[iObj]) ;
3694                    if (objSOS)
3695                      continue;
3696                    int iColumn = objects[iObj]->columnNumber();
3697                    assert (iColumn>=0);
3698                    if (originalColumns)
3699                      iColumn = originalColumns[iColumn];
3700                    if (branchDirection) {
3701                      CbcSimpleInteger * obj =
3702                        dynamic_cast <CbcSimpleInteger *>(objects[iObj]) ;
3703                      if (obj) { 
3704                        obj->setPreferredWay(branchDirection[iColumn]);
3705                      } else {
3706                        CbcObject * obj =
3707                          dynamic_cast <CbcObject *>(objects[iObj]) ;
3708                        assert (obj);
3709                        obj->setPreferredWay(branchDirection[iColumn]);
3710                      }
3711                    }
3712                    if (priorities) {
3713                      int iPriority = priorities[iColumn];
3714                      if (iPriority>0)
3715                        objects[iObj]->setPriority(iPriority);
3716                    }
3717                    if (pseudoUp&&pseudoUp[iColumn]) {
3718                      CbcSimpleIntegerPseudoCost * obj1a =
3719                        dynamic_cast <CbcSimpleIntegerPseudoCost *>(objects[iObj]) ;
3720                      assert (obj1a);
3721                      if (pseudoDown[iColumn]>0.0)
3722                        obj1a->setDownPseudoCost(pseudoDown[iColumn]);
3723                      if (pseudoUp[iColumn]>0.0)
3724                        obj1a->setUpPseudoCost(pseudoUp[iColumn]);
3725                    }
3726                  }
3727                  // *************************************************************
3728                } else {
3729                  // *************************************************************
3730                  // OsiObjects
3731                  // Find if none
3732                  int numberIntegers = testOsiSolver->getNumIntegers();
3733                  /* model may not have created objects
3734                     If none then create
3735                  */
3736                  if (!numberIntegers||!testOsiSolver->numberObjects()) {
3737                    //int type = (pseudoUp) ? 1 : 0;
3738                    testOsiSolver->findIntegers(false);
3739                    numberIntegers = testOsiSolver->getNumIntegers();
3740                  }
3741                  if (preProcess&&process.numberSOS()) {
3742                    int numberSOS = process.numberSOS();
3743                    OsiObject ** oldObjects = testOsiSolver->objects();
3744                    // Do sets and priorities
3745                    OsiObject ** objects = new OsiObject * [numberSOS];
3746                    // set old objects to have low priority
3747                    int numberOldObjects = testOsiSolver->numberObjects();
3748                    int numberColumns = testOsiSolver->getNumCols();
3749                    for (int iObj = 0;iObj<numberOldObjects;iObj++) {
3750                      oldObjects[iObj]->setPriority(numberColumns+1);
3751                      int iColumn = oldObjects[iObj]->columnNumber();
3752                      assert (iColumn>=0);
3753                      if (iColumn>=numberOriginalColumns)
3754                        continue;
3755                      if (originalColumns)
3756                        iColumn = originalColumns[iColumn];
3757                      if (branchDirection) {
3758                        OsiSimpleInteger * obj =
3759                          dynamic_cast <OsiSimpleInteger *>(oldObjects[iObj]) ;
3760                        if (obj) { 
3761                          obj->setPreferredWay(branchDirection[iColumn]);
3762                        } else {
3763                          OsiObject2 * obj =
3764                            dynamic_cast <OsiObject2 *>(oldObjects[iObj]) ;
3765                          if (obj)
3766                            obj->setPreferredWay(branchDirection[iColumn]);
3767                        }
3768                      }
3769                      if (pseudoUp) {
3770                        abort();
3771                      }
3772                    }
3773                    const int * starts = process.startSOS();
3774                    const int * which = process.whichSOS();
3775                    const int * type = process.typeSOS();
3776                    const double * weight = process.weightSOS();
3777                    int iSOS;
3778                    for (iSOS =0;iSOS<numberSOS;iSOS++) {
3779                      int iStart = starts[iSOS];
3780                      int n=starts[iSOS+1]-iStart;
3781                      objects[iSOS] = new OsiSOS(testOsiSolver,n,which+iStart,weight+iStart,
3782                                                 type[iSOS]);
3783                      // branch on long sets first
3784                      objects[iSOS]->setPriority(numberColumns-n);
3785                    }
3786                    testOsiSolver->addObjects(numberSOS,objects);
3787                    for (iSOS=0;iSOS<numberSOS;iSOS++)
3788                      delete objects[iSOS];
3789                    delete [] objects;
3790                  } else if (priorities||branchDirection||pseudoDown||pseudoUp||numberSOS) {
3791                    if (numberSOS) {
3792                      // Do sets and priorities
3793                      OsiObject ** objects = new OsiObject * [numberSOS];
3794                      int iSOS;
3795                      if (originalColumns) {
3796                        // redo sequence numbers
3797                        int numberColumns = testOsiSolver->getNumCols();
3798                        int nOld = originalColumns[numberColumns-1]+1;
3799                        int * back = new int[nOld];
3800                        int i;
3801                        for (i=0;i<nOld;i++)
3802                          back[i]=-1;
3803                        for (i=0;i<numberColumns;i++)
3804                          back[originalColumns[i]]=i;
3805                        // Really need better checks
3806                        int nMissing=0;
3807                        int n=sosStart[numberSOS];
3808                        for (i=0;i<n;i++) {
3809                          int iColumn = sosIndices[i];
3810                          int jColumn = back[iColumn];
3811                          if (jColumn>=0) 
3812                            sosIndices[i] = jColumn;
3813                          else 
3814                            nMissing++;
3815                        }
3816                        delete [] back;
3817                        if (nMissing)
3818                          printf("%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
3819                      }
3820                      for (iSOS =0;iSOS<numberSOS;iSOS++) {
3821                        int iStart = sosStart[iSOS];
3822                        int n=sosStart[iSOS+1]-iStart;
3823                        objects[iSOS] = new OsiSOS(testOsiSolver,n,sosIndices+iStart,sosReference+iStart,
3824                                                   sosType[iSOS]);
3825                        if (sosPriority)
3826                          objects[iSOS]->setPriority(sosPriority[iSOS]);
3827                        else if (!prioritiesIn)
3828                          objects[iSOS]->setPriority(10);  // rather than 1000
3829                      }
3830                      // delete any existing SOS objects
3831                      int numberObjects=testOsiSolver->numberObjects();
3832                      OsiObject ** oldObjects=testOsiSolver->objects();
3833                      int nNew=0;
3834                      for (int i=0;i<numberObjects;i++) {
3835                        OsiObject * objThis = oldObjects[i];
3836                        OsiSOS * obj1 =
3837                          dynamic_cast <OsiSOS *>(objThis) ;
3838                        OsiSOS * obj2 =
3839                          dynamic_cast <OsiSOS *>(objThis) ;
3840                        if (!obj1&&!obj2) {
3841                          oldObjects[nNew++]=objThis;
3842                        } else {
3843                          delete objThis;
3844                        }
3845                      }
3846                      testOsiSolver->setNumberObjects(nNew);
3847                      testOsiSolver->addObjects(numberSOS,objects);
3848                      for (iSOS=0;iSOS<numberSOS;iSOS++)
3849                        delete objects[iSOS];
3850                      delete [] objects;
3851                    }
3852                  }
3853                  OsiObject ** objects = testOsiSolver->objects();
3854                  int numberObjects = testOsiSolver->numberObjects();
3855                  int logLevel = parameters[log].intValue();
3856                  for (int iObj = 0;iObj<numberObjects;iObj++) {
3857                    // skip sos
3858                    OsiSOS * objSOS =
3859                      dynamic_cast <OsiSOS *>(objects[iObj]) ;
3860                    if (objSOS) {
3861                      if (logLevel>2)
3862                        printf("Set %d is SOS - priority %d\n",iObj,objSOS->priority());
3863                      continue;
3864                    }
3865                    int iColumn = objects[iObj]->columnNumber();
3866                    if (iColumn>=0) {
3867                      if (originalColumns)
3868                        iColumn = originalColumns[iColumn];
3869                      if (branchDirection) {
3870                        OsiSimpleInteger * obj =
3871                          dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
3872                        if (obj) { 
3873                          obj->setPreferredWay(branchDirection[iColumn]);
3874                        } else {
3875                          OsiObject2 * obj =
3876                            dynamic_cast <OsiObject2 *>(objects[iObj]) ;
3877                          if (obj)
3878                            obj->setPreferredWay(branchDirection[iColumn]);
3879                        }
3880                      }
3881                    }
3882                    if (priorities) {
3883                      int iPriority = priorities[iColumn];
3884                      if (iPriority>0)
3885                        objects[iObj]->setPriority(iPriority);
3886                    }
3887                    if (logLevel>2)
3888                      printf("Obj %d is int? - priority %d\n",iObj,objects[iObj]->priority());
3889                    if (pseudoUp&&pseudoUp[iColumn]) {
3890                      abort();
3891                    }
3892                  }
3893                  // *************************************************************
3894                }
3895                int statistics = (printOptions>0) ? printOptions: 0;
3896#ifdef COIN_HAS_ASL
3897                if (!usingAmpl) {
3898#endif
3899                  free(priorities);
3900                  priorities=NULL;
3901                  free(branchDirection);
3902                  branchDirection=NULL;
3903                  free(pseudoDown);
3904                  pseudoDown=NULL;
3905                  free(pseudoUp);
3906                  pseudoUp=NULL;
3907                  free(solutionIn);
3908                  solutionIn=NULL;
3909                  free(prioritiesIn);
3910                  prioritiesIn=NULL;
3911                  free(sosStart);
3912                  sosStart=NULL;
3913                  free(sosIndices);
3914                  sosIndices=NULL;
3915                  free(sosType);
3916                  sosType=NULL;
3917                  free(sosReference);
3918                  sosReference=NULL;
3919                  free(cut);
3920                  cut=NULL;
3921                  free(sosPriority);
3922                  sosPriority=NULL;
3923#ifdef COIN_HAS_ASL
3924                }
3925#endif               
3926                if (nodeStrategy) {
3927                  // change default
3928                  if (nodeStrategy>1) {
3929                    // up or down
3930                    int way = ((nodeStrategy%1)==1) ? -1 : +1;
3931                    babModel->setPreferredWay(way);
3932#if 0
3933                    OsiObject ** objects = babModel->objects();
3934                    int numberObjects = babModel->numberObjects();
3935                    for (int iObj = 0;iObj<numberObjects;iObj++) {
3936                      CbcObject * obj =
3937                        dynamic_cast <CbcObject *>(objects[iObj]) ;
3938                      assert (obj);
3939                      obj->setPreferredWay(way);
3940                    }
3941#endif
3942                  }
3943                  if (nodeStrategy==1||nodeStrategy>3) {
3944                    // depth
3945                    CbcCompareDefault compare;
3946                    compare.setWeight(-3.0);
3947                    babModel->setNodeComparison(compare);
3948                  }
3949                }
3950                if (cppValue>=0) {
3951                  int prepro = useStrategy ? -1 : preProcess;
3952                  // generate code
3953                  FILE * fp = fopen("user_driver.cpp","w");
3954                  if (fp) {
3955                    // generate enough to do BAB
3956                    babModel->generateCpp(fp,1);
3957                    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
3958                    // Make general so do factorization
3959                    int factor = osiclp->getModelPtr()->factorizationFrequency();
3960                    osiclp->getModelPtr()->setFactorizationFrequency(200);
3961                    osiclp->generateCpp(fp);
3962                    osiclp->getModelPtr()->setFactorizationFrequency(factor);
3963                    //solveOptions.generateCpp(fp);
3964                    fclose(fp);
3965                    // now call generate code
3966                    generateCode(babModel,"user_driver.cpp",cppValue,prepro);
3967                  } else {
3968                    std::cout<<"Unable to open file user_driver.cpp"<<std::endl;
3969                  }
3970                }
3971                if (!babModel->numberStrong())
3972                  babModel->setNumberBeforeTrust(0);
3973                if (useStrategy) {
3974                  CbcStrategyDefault strategy(true,babModel->numberStrong(),babModel->numberBeforeTrust());
3975                  strategy.setupPreProcessing(1);
3976                  babModel->setStrategy(strategy);
3977                }
3978                if (testOsiOptions>=0) {
3979                  printf("Testing OsiObject options %d\n",testOsiOptions);
3980                  if (!numberSOS) {
3981                    babModel->solver()->findIntegersAndSOS(false);
3982#ifdef COIN_HAS_LINK
3983                    // If linked then pass in model
3984                    OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
3985                    if (solver3) {
3986                      int options = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()/10000;
3987                      CglStored stored;
3988                      if (options) {
3989                        printf("nlp options %d\n",options);
3990                        /*
3991                          1 - force mini branch and bound
3992                          2 - set priorities high on continuous
3993                          4 - try adding OA cuts
3994                          8 - try doing quadratic linearization
3995                          16 - try expanding knapsacks
3996                          32 - OA cuts strictly concave
3997                        */
3998                        if ((options&2)) {
3999                          solver3->setBiLinearPriorities(10,tightenFactor > 0.0 ? tightenFactor : 1.0);
4000                        } else if (tightenFactor>0.0) {
4001                          // set grid size for all continuous bi-linear
4002                          solver3->setMeshSizes(tightenFactor);
4003                        }
4004                        if ((options&4)) {
4005                          solver3->setSpecialOptions2(solver3->specialOptions2()|(8+4));
4006                          // say convex
4007                          solver3->sayConvex((options&32)==0);
4008                        }
4009                        if ((options&1)!=0&&slpValue>0)
4010                          solver3->setFixedPriority(slpValue);
4011                        double cutoff=COIN_DBL_MAX;
4012                        if ((options&8))
4013                          cutoff=solver3->linearizedBAB(&stored);
4014                        if (cutoff<babModel->getCutoff())
4015                          babModel->setCutoff(cutoff);
4016                      }
4017                      solver3->setCbcModel(babModel);
4018                      babModel->addCutGenerator(&stored,1,"Stored");
4019                      CglTemporary temp;
4020                      babModel->addCutGenerator(&temp,1,"OnceOnly");
4021                      //choose.setNumberBeforeTrusted(2000);
4022                      //choose.setNumberStrong(20);
4023                    }
4024                    // For temporary testing of heuristics
4025                    //int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
4026                    if (testOsiOptions>=10) {
4027                      printf("*** Temp heuristic with mode %d\n",testOsiOptions-10);
4028                      OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
4029                      assert (solver3) ;
4030                      double * solution = solver3->heuristicSolution(slpValue>0 ? slpValue : 20 ,1.0e-5,testOsiOptions-10);
4031                      assert(solution);
4032                    }
4033#endif
4034                  } else {
4035                    // move across
4036                    babModel->deleteObjects(false);
4037                    //babModel->addObjects(babModel->solver()->numberObjects(),babModel->solver()->objects());
4038                  }
4039                  CbcBranchDefaultDecision decision;
4040                  if (babModel->numberStrong()) {
4041                    OsiChooseStrong choose(babModel->solver());
4042                    choose.setNumberBeforeTrusted(babModel->numberBeforeTrust());
4043                    choose.setNumberStrong(babModel->numberStrong());
4044                    choose.setShadowPriceMode(testOsiOptions);
4045                    decision.setChooseMethod(choose);
4046                  } else {
4047                    OsiChooseVariable choose(babModel->solver());
4048                    decision.setChooseMethod(choose);
4049                  }
4050                  babModel->setBranchingMethod(decision);
4051                  if (useCosts&&testOsiOptions>=0) {
4052                    int numberColumns = babModel->getNumCols();
4053                    int * sort = new int[numberColumns];
4054                    double * dsort = new double[numberColumns];
4055                    int * priority = new int [numberColumns];
4056                    const double * objective = babModel->getObjCoefficients();
4057                    const double * lower = babModel->getColLower() ;
4058                    const double * upper = babModel->getColUpper() ;
4059                    const CoinPackedMatrix * matrix = babModel->solver()->getMatrixByCol();
4060                    const int * columnLength = matrix->getVectorLengths();
4061                    int iColumn;
4062                    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4063                      sort[iColumn]=iColumn;
4064                      if (useCosts==1)
4065                        dsort[iColumn]=-fabs(objective[iColumn]);
4066                      else if (useCosts==2)
4067                        dsort[iColumn]=iColumn;
4068                      else if (useCosts==3)
4069                        dsort[iColumn]=upper[iColumn]-lower[iColumn];
4070                      else if (useCosts==4)
4071                        dsort[iColumn]=-(upper[iColumn]-lower[iColumn]);
4072                      else if (useCosts==5)
4073                        dsort[iColumn]=-columnLength[iColumn];
4074                    }
4075                    CoinSort_2(dsort,dsort+numberColumns,sort);
4076                    int level=0;
4077                    double last = -1.0e100;
4078                    for (int i=0;i<numberColumns;i++) {
4079                      int iPut=sort[i];
4080                      if (dsort[i]!=last) {
4081                        level++;
4082                        last=dsort[i];
4083                      }
4084                      priority[iPut]=level;
4085                    }
4086                    OsiObject ** objects = babModel->objects();
4087                    int numberObjects = babModel->numberObjects();
4088                    for (int iObj = 0;iObj<numberObjects;iObj++) {
4089                      OsiObject * obj = objects[iObj] ;
4090                      int iColumn = obj->columnNumber();
4091                      if (iColumn>=0)
4092                        obj->setPriority(priority[iColumn]);
4093                    }
4094                    delete [] priority;
4095                    delete [] sort;
4096                    delete [] dsort;
4097                  }
4098                }
4099                checkSOS(babModel, babModel->solver());
4100                if (doSprint>0) {
4101                  // Sprint for primal solves
4102                  ClpSolve::SolveType method = ClpSolve::usePrimalorSprint;
4103                  ClpSolve::PresolveType presolveType = ClpSolve::presolveOff;
4104                  int numberPasses = 5;
4105                  int options[] = {0,3,0,0,0,0};
4106                  int extraInfo[] = {-1,20,-1,-1,-1,-1};
4107                  extraInfo[1]=doSprint;
4108                  int independentOptions[] = {0,0,3};
4109                  ClpSolve clpSolve(method,presolveType,numberPasses,
4110                                    options,extraInfo,independentOptions);
4111                  // say use in OsiClp
4112                  clpSolve.setSpecialOption(6,1);
4113                  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4114                  osiclp->setSolveOptions(clpSolve);
4115                  osiclp->setHintParam(OsiDoDualInResolve,false);
4116                  // switch off row copy
4117                  osiclp->getModelPtr()->setSpecialOptions(osiclp->getModelPtr()->specialOptions()|256);
4118                  osiclp->getModelPtr()->setInfeasibilityCost(1.0e11);
4119                }
4120#ifdef COIN_HAS_LINK
4121                if (storedAmpl.sizeRowCuts()) {
4122                  if (preProcess) {
4123                    const int * originalColumns = process.originalColumns();
4124                    int numberColumns = babModel->getNumCols();
4125                    int * newColumn = new int[numberOriginalColumns];
4126                    int i;
4127                    for (i=0;i<numberOriginalColumns;i++) 
4128                      newColumn[i]=-1;
4129                    for (i=0;i<numberColumns;i++) {
4130                      int iColumn = originalColumns[i];
4131                      newColumn[iColumn]=i;
4132                    }
4133                    int * buildColumn = new int[numberColumns];
4134                    // Build up valid cuts
4135                    int nBad=0;
4136                    int nCuts = storedAmpl.sizeRowCuts();
4137                    CglStored newCuts;
4138                    for (i=0;i<nCuts;i++) {
4139                      const OsiRowCut * cut = storedAmpl.rowCutPointer(i);
4140                      double lb = cut->lb();
4141                      double ub = cut->ub();
4142                      int n=cut->row().getNumElements();
4143                      const int * column = cut->row().getIndices();
4144                      const double * element = cut->row().getElements();
4145                      bool bad=false;
4146                      for (int i=0;i<n;i++) {
4147                        int iColumn = column[i];
4148                        iColumn = newColumn[iColumn];
4149                        if (iColumn>=0) {
4150                          buildColumn[i]=iColumn;
4151                        } else {
4152                          bad=true;
4153                          break;
4154                        }
4155                      }
4156                      if (!bad) {
4157                        newCuts.addCut(lb,ub,n,buildColumn,element);
4158                      } else {
4159                        nBad++;
4160                      }
4161                    }
4162                    storedAmpl=newCuts;
4163                    if (nBad)
4164                      printf("%d cuts dropped\n",nBad);
4165                    delete [] newColumn;
4166                    delete [] buildColumn;
4167                  }
4168                }
4169#endif
4170#ifdef CLP_MALLOC_STATISTICS
4171                malloc_stats();
4172                malloc_stats2();
4173#endif
4174                if (outputFormat==5) {
4175                  osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4176                  lpSolver = osiclp->getModelPtr();
4177                  lpSolver->setPersistenceFlag(1);
4178                }
4179#ifdef COIN_HAS_ASL
4180                // add in lotsizing
4181                if (usingAmpl&&info.special) {
4182                  int numberColumns = babModel->getNumCols();
4183                  int i;
4184                  int n=0;
4185                  if (preProcess) {
4186                    const int * originalColumns = process.originalColumns();
4187                    for (i=0;i<numberColumns;i++) {
4188                      int iColumn = originalColumns[i];
4189                      assert (iColumn>=i);
4190                      int iType = info.special[iColumn];
4191                      if (iType) {
4192                        assert (iType==1);
4193                        n++;
4194                      }
4195                      info.special[i]=iType;
4196                    }
4197                  }
4198                  if (n) {
4199                    int numberIntegers=0;
4200                    int numberOldObjects=0;
4201                    OsiObject ** oldObjects=NULL;
4202                    const double * lower = babModel->solver()->getColLower();
4203                    const double * upper = babModel->solver()->getColUpper();
4204                    if (testOsiOptions<0) {
4205                      // *************************************************************
4206                      // CbcObjects
4207                      numberIntegers = babModel->numberIntegers();
4208                      /* model may not have created objects
4209                         If none then create
4210                      */
4211                      if (!numberIntegers||!babModel->numberObjects()) {
4212                        int type = (pseudoUp) ? 1 : 0;
4213                        babModel->findIntegers(true,type);
4214                        numberIntegers = babModel->numberIntegers();
4215                      }
4216                      oldObjects = babModel->objects();
4217                      numberOldObjects = babModel->numberObjects();
4218                    } else {
4219                      numberIntegers = testOsiSolver->getNumIntegers();
4220                      if (!numberIntegers||!testOsiSolver->numberObjects()) {
4221                        /* model may not have created objects
4222                           If none then create
4223                        */
4224                        testOsiSolver->findIntegers(false);
4225                        numberIntegers = testOsiSolver->getNumIntegers();
4226                      }
4227                      oldObjects = testOsiSolver->objects();
4228                      numberOldObjects = testOsiSolver->numberObjects();
4229                    }
4230                    OsiObject ** objects = new OsiObject * [n];
4231                    n=0;
4232                    // set new objects to have one lower priority
4233                    double ranges[] = {-COIN_DBL_MAX,-1.0,1.0,COIN_DBL_MAX};
4234                    for (int iObj = 0;iObj<numberOldObjects;iObj++) {
4235                      int iColumn = oldObjects[iObj]->columnNumber();
4236                      if (iColumn>=0&&info.special[iColumn]) {
4237                        if (lower[iColumn]<=-1.0&&upper[iColumn]>=0.0) {
4238                          ranges[0]=lower[iColumn];
4239                          ranges[3]=upper[iColumn];
4240                          int priority = oldObjects[iObj]->priority();
4241                          if (testOsiOptions<0) {
4242                            objects[n] = new CbcLotsize(babModel,iColumn,2,ranges,true);
4243                          } else {
4244                            objects[n] = new OsiLotsize(testOsiSolver,iColumn,2,ranges,true);
4245                          }
4246                          objects[n++]->setPriority (priority-1);
4247                        }
4248                      }
4249                    }
4250                    if (testOsiOptions<0) {
4251                      babModel->addObjects(n,objects);
4252                    } else {
4253                      testOsiSolver->addObjects(n,objects);
4254                    }
4255                    for (i=0;i<n;i++)
4256                      delete objects[i];
4257                    delete [] objects;
4258                  }
4259                }
4260                if (storedAmpl.sizeRowCuts()) {
4261                  //babModel->addCutGenerator(&storedAmpl,1,"AmplStored");
4262                  int numberRowCuts = storedAmpl.sizeRowCuts();
4263                  for (int i=0;i<numberRowCuts;i++) {
4264                    const OsiRowCut * rowCutPointer = storedAmpl.rowCutPointer(i);
4265                    babModel->makeGlobalCut(rowCutPointer);
4266                  }
4267                }
4268#endif
4269                babModel->branchAndBound(statistics);
4270#ifdef CLP_MALLOC_STATISTICS
4271                malloc_stats();
4272                malloc_stats2();
4273#endif
4274                checkSOS(babModel, babModel->solver());
4275              } else if (type==MIPLIB) {
4276                CbcStrategyDefault strategy(true,babModel->numberStrong(),babModel->numberBeforeTrust());
4277                // Set up pre-processing
4278                int translate2[]={9999,1,1,3,2,4,5};
4279                if (preProcess)
4280                  strategy.setupPreProcessing(translate2[preProcess]);
4281                babModel->setStrategy(strategy);
4282                if (outputFormat==5) {
4283                  osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4284                  lpSolver = osiclp->getModelPtr();
4285                  lpSolver->setPersistenceFlag(1);
4286                }
4287                if (testOsiOptions>=0) {
4288                  printf("Testing OsiObject options %d\n",testOsiOptions);
4289                  CbcBranchDefaultDecision decision;
4290                  OsiChooseStrong choose(babModel->solver());
4291                  choose.setNumberBeforeTrusted(babModel->numberBeforeTrust());
4292                  choose.setNumberStrong(babModel->numberStrong());
4293                  choose.setShadowPriceMode(testOsiOptions);
4294                  //babModel->deleteObjects(false);
4295                  decision.setChooseMethod(choose);
4296                  babModel->setBranchingMethod(decision);
4297                }
4298                *miplibSolver = babModel;
4299                return 777;
4300              } else {
4301                strengthenedModel = babModel->strengthenedModel();
4302              }
4303              currentBranchModel = NULL;
4304              osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4305              if (debugFile=="createAfterPre"&&babModel->bestSolution()) {
4306                lpSolver = osiclp->getModelPtr();
4307                //move best solution (should be there -- but ..)
4308                int n = lpSolver->getNumCols();
4309                memcpy(lpSolver->primalColumnSolution(),babModel->bestSolution(),n*sizeof(double));
4310                saveSolution(osiclp->getModelPtr(),"debug.file");
4311              }
4312              if (!noPrinting) {
4313                // Print more statistics
4314                std::cout<<"Cuts at root node changed objective from "<<babModel->getContinuousObjective()
4315                         <<" to "<<babModel->rootObjectiveAfterCuts()<<std::endl;
4316               
4317                numberGenerators = babModel->numberCutGenerators();
4318                for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
4319                  CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
4320                  std::cout<<generator->cutGeneratorName()<<" was tried "
4321                           <<generator->numberTimesEntered()<<" times and created "
4322                           <<generator->numberCutsInTotal()<<" cuts of which "
4323                           <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
4324                  if (generator->timing())
4325                    std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
4326                  else
4327                    std::cout<<std::endl;
4328                }
4329              }
4330              time2 = CoinCpuTime();
4331              totalTime += time2-time1;
4332              // For best solution
4333              double * bestSolution = NULL;
4334              if (babModel->getMinimizationObjValue()<1.0e50&&type==BAB) {
4335                // post process
4336                if (preProcess) {
4337                  int n = saveSolver->getNumCols();
4338                  bestSolution = new double [n];
4339                  process.postProcess(*babModel->solver());
4340                  // Solution now back in saveSolver
4341                  babModel->assignSolver(saveSolver);
4342                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
4343                } else {
4344                  int n = babModel->solver()->getNumCols();
4345                  bestSolution = new double [n];
4346                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
4347                }
4348                // and put back in very original solver
4349                {
4350                  ClpSimplex * original = originalSolver.getModelPtr();
4351                  double * lower = original->columnLower();
4352                  double * upper = original->columnUpper();
4353                  double * solution = original->primalColumnSolution();
4354                  int n = original->numberColumns();
4355                  assert (!n||n==babModel->solver()->getNumCols());
4356                  for (int i=0;i<n;i++) {
4357                    solution[i]=bestSolution[i];
4358                    if (solver1.isInteger(i)) {
4359                      lower[i]=solution[i];
4360                      upper[i]=solution[i];
4361                    }
4362                  }
4363                }
4364                checkSOS(babModel, babModel->solver());
4365              }
4366              if (type==STRENGTHEN&&strengthenedModel)
4367                clpSolver = dynamic_cast< OsiClpSolverInterface*> (strengthenedModel);
4368#ifdef COIN_HAS_ASL
4369              else if (usingAmpl) 
4370                clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4371#endif
4372              lpSolver = clpSolver->getModelPtr();
4373              if (numberChanged) {
4374                for (int i=0;i<numberChanged;i++) {
4375                  int iColumn=changed[i];
4376                  clpSolver->setContinuous(iColumn);
4377                }
4378                delete [] changed;
4379              }
4380              if (type==BAB) {
4381                //move best solution (should be there -- but ..)
4382                int n = lpSolver->getNumCols();
4383                if (bestSolution) {
4384                  memcpy(lpSolver->primalColumnSolution(),bestSolution,n*sizeof(double));
4385                  // now see what that does to row solution
4386                  int numberRows=lpSolver->numberRows();
4387                  double * rowSolution = lpSolver->primalRowSolution();
4388                  memset (rowSolution,0,numberRows*sizeof(double));
4389                  lpSolver->clpMatrix()->times(1.0,bestSolution,rowSolution);
4390                }
4391                if (debugFile=="create"&&bestSolution) {
4392                  saveSolution(lpSolver,"debug.file");
4393                }
4394                delete [] bestSolution;
4395                std::string statusName[]={"Finished","Stopped on ","Difficulties",
4396                                          "","","User ctrl-c"};
4397                std::string minor[]={"","","gap","nodes","time","","solutions","user ctrl-c"};
4398                int iStat = babModel->status();
4399                int iStat2 = babModel->secondaryStatus();
4400                if (!noPrinting)
4401                  std::cout<<"Result - "<<statusName[iStat]<<minor[iStat2]
4402                           <<" objective "<<babModel->getObjValue()<<
4403                    " after "<<babModel->getNodeCount()<<" nodes and "
4404                           <<babModel->getIterationCount()<<
4405                    " iterations - took "<<time2-time1<<" seconds"<<std::endl;
4406#ifdef COIN_HAS_ASL
4407                if (usingAmpl) {
4408                  clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4409                  lpSolver = clpSolver->getModelPtr();
4410                  double value = babModel->getObjValue()*lpSolver->getObjSense();
4411                  char buf[300];
4412                  int pos=0;
4413                  if (iStat==0) {
4414                    if (babModel->getObjValue()<1.0e40) {
4415                      pos += sprintf(buf+pos,"optimal," );
4416                    } else {
4417                      // infeasible
4418                      iStat=1;
4419                      pos += sprintf(buf+pos,"infeasible,");
4420                    }
4421                  } else if (iStat==1) {
4422                    if (iStat2!=6)
4423                      iStat=3;
4424                    else
4425                      iStat=4;
4426                    pos += sprintf(buf+pos,"stopped on %s,",minor[iStat2].c_str());
4427                  } else if (iStat==2) {
4428                    iStat = 7;
4429                    pos += sprintf(buf+pos,"stopped on difficulties,");
4430                  } else if (iStat==5) {
4431                    iStat = 3;
4432                    pos += sprintf(buf+pos,"stopped on ctrl-c,");
4433                  } else {
4434                    pos += sprintf(buf+pos,"status unknown,");
4435                    iStat=6;
4436                  }
4437                  info.problemStatus=iStat;
4438                  info.objValue = value;
4439                  if (babModel->getObjValue()<1.0e40) {
4440                    int precision = ampl_obj_prec();
4441                    if (precision>0)
4442                      pos += sprintf(buf+pos," objective %.*g",precision,
4443                                     value);
4444                    else
4445                      pos += sprintf(buf+pos," objective %g",value);
4446                  }
4447                  sprintf(buf+pos,"\n%d nodes, %d iterations, %g seconds",
4448                          babModel->getNodeCount(),
4449                          babModel->getIterationCount(),
4450                          totalTime);
4451                  if (bestSolution) {
4452                    free(info.primalSolution);
4453                    if (!numberKnapsack) {
4454                      info.primalSolution = (double *) malloc(n*sizeof(double));
4455                      CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution);
4456                      int numberRows = lpSolver->numberRows();
4457                      free(info.dualSolution);
4458                      info.dualSolution = (double *) malloc(numberRows*sizeof(double));
4459                      CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution);
4460                    } else {
4461                      // expanded knapsack
4462                      info.dualSolution=NULL;
4463                      int numberColumns = saveCoinModel.numberColumns();
4464                      info.primalSolution = (double *) malloc(numberColumns*sizeof(double));
4465                      // Fills in original solution (coinModel length)
4466                      afterKnapsack(saveTightenedModel,  whichColumn,  knapsackStart, 
4467                                    knapsackRow,  numberKnapsack,
4468                                    lpSolver->primalColumnSolution(), info.primalSolution,1);
4469                    }
4470                  } else {
4471                    info.primalSolution=NULL;
4472                    info.dualSolution=NULL;
4473                  }
4474                  // put buffer into info
4475                  strcpy(info.buffer,buf);
4476                }
4477#endif
4478              } else {
4479                std::cout<<"Model strengthened - now has "<<clpSolver->getNumRows()
4480                         <<" rows"<<std::endl;
4481              }
4482              time1 = time2;
4483              delete babModel;
4484              babModel=NULL;
4485            } else {
4486              std::cout << "** Current model not valid" << std::endl ; 
4487            }
4488            break ;
4489          case IMPORT:
4490            {
4491#ifdef COIN_HAS_ASL
4492              if (!usingAmpl) {
4493#endif
4494                free(priorities);
4495                priorities=NULL;
4496                free(branchDirection);
4497                branchDirection=NULL;
4498                free(pseudoDown);
4499                pseudoDown=NULL;
4500                free(pseudoUp);
4501                pseudoUp=NULL;
4502                free(solutionIn);
4503                solutionIn=NULL;
4504                free(prioritiesIn);
4505                prioritiesIn=NULL;
4506                free(sosStart);
4507                sosStart=NULL;
4508                free(sosIndices);
4509                sosIndices=NULL;
4510                free(sosType);
4511                sosType=NULL;
4512                free(sosReference);
4513                sosReference=NULL;
4514                free(cut);
4515                cut=NULL;
4516                free(sosPriority);
4517                sosPriority=NULL;
4518#ifdef COIN_HAS_ASL
4519              }
4520#endif               
4521              delete babModel;
4522              babModel=NULL;
4523              // get next field
4524              field = CoinReadGetString(argc,argv);
4525              if (field=="$") {
4526                field = parameters[iParam].stringValue();
4527              } else if (field=="EOL") {
4528                parameters[iParam].printString();
4529                break;
4530              } else {
4531                parameters[iParam].setStringValue(field);
4532              }
4533              std::string fileName;
4534              bool canOpen=false;
4535              // See if gmpl file
4536              int gmpl=0;
4537              std::string gmplData;
4538              if (field=="-") {
4539                // stdin
4540                canOpen=true;
4541                fileName = "-";
4542              } else {
4543                // See if .lp
4544                {
4545                  const char * c_name = field.c_str();
4546                  int length = strlen(c_name);
4547                  if (length>3&&!strncmp(c_name+length-3,".lp",3))
4548                    gmpl=-1; // .lp
4549                }
4550                bool absolutePath;
4551                if (dirsep=='/') {
4552                  // non Windows (or cygwin)
4553                  absolutePath=(field[0]=='/');
4554                } else {
4555                  //Windows (non cycgwin)
4556                  absolutePath=(field[0]=='\\');
4557                  // but allow for :
4558                  if (strchr(field.c_str(),':'))
4559                    absolutePath=true;
4560                }
4561                if (absolutePath) {
4562                  fileName = field;
4563                } else if (field[0]=='~') {
4564                  char * environVar = getenv("HOME");
4565                  if (environVar) {
4566                    std::string home(environVar);
4567                    field=field.erase(0,1);
4568                    fileName = home+field;
4569                  } else {
4570                    fileName=field;
4571                  }
4572                } else {
4573                  fileName = directory+field;
4574                  // See if gmpl (model & data) - or even lp file
4575                  int length = field.size();
4576                  int percent = field.find('%');
4577                  if (percent<length&&percent>0) {
4578                    gmpl=1;
4579                    fileName = directory+field.substr(0,percent);
4580                    gmplData = directory+field.substr(percent+1);
4581                    if (percent<length-1)
4582                      gmpl=2; // two files
4583                    printf("GMPL model file %s and data file %s\n",
4584                           fileName.c_str(),gmplData.c_str());
4585                  }
4586                }
4587                std::string name=fileName;
4588                if (fileCoinReadable(name)) {
4589                  // can open - lets go for it
4590                  canOpen=true;
4591                  if (gmpl==2) {
4592                    FILE *fp;
4593                    fp=fopen(gmplData.c_str(),"r");
4594                    if (fp) {
4595                      fclose(fp);
4596                    } else {
4597                      canOpen=false;
4598                      std::cout<<"Unable to open file "<<gmplData<<std::endl;
4599                    }
4600                  }
4601                } else {
4602                  std::cout<<"Unable to open file "<<fileName<<std::endl;
4603                }
4604              }
4605              if (canOpen) {
4606                int status;
4607                ClpSimplex * lpSolver = clpSolver->getModelPtr();
4608                if (!gmpl)
4609                  status =clpSolver->readMps(fileName.c_str(),
4610                                                 keepImportNames!=0,
4611                                                 allowImportErrors!=0);
4612                else if (gmpl>0)
4613                  status= lpSolver->readGMPL(fileName.c_str(),
4614                                                  (gmpl==2) ? gmplData.c_str() : NULL,
4615                                                  keepImportNames!=0);
4616                else
4617                  status= lpSolver->readLp(fileName.c_str(),1.0e-12);
4618                if (!status||(status>0&&allowImportErrors)) {
4619                  goodModel=true;
4620                  // sets to all slack (not necessary?)
4621                  lpSolver->createStatus();
4622                  // make sure integer
4623                  int numberColumns = lpSolver->numberColumns();
4624                  for (int i=0;i<numberColumns;i++) {
4625                    if (lpSolver->isInteger(i))
4626                      clpSolver->setInteger(i);
4627                  }
4628                  time2 = CoinCpuTime();
4629                  totalTime += time2-time1;
4630                  time1=time2;
4631                  // Go to canned file if just input file
4632                  if (CbcOrClpRead_mode==2&&argc==2) {
4633                    // only if ends .mps
4634                    char * find = (char *)strstr(fileName.c_str(),".mps");
4635                    if (find&&find[4]=='\0') {
4636                      find[1]='p'; find[2]='a';find[3]='r';
4637                      FILE *fp=fopen(fileName.c_str(),"r");
4638                      if (fp) {
4639                        CbcOrClpReadCommand=fp; // Read from that file
4640                        CbcOrClpRead_mode=-1;
4641                      }
4642                    }
4643                  }
4644                } else {
4645                  // errors
4646                  std::cout<<"There were "<<status<<
4647                    " errors on input"<<std::endl;
4648                }
4649              }
4650            }
4651            break;
4652          case MODELIN:
4653#ifdef COIN_HAS_LINK
4654            {
4655              // get next field
4656              field = CoinReadGetString(argc,argv);
4657              if (field=="$") {
4658                field = parameters[iParam].stringValue();
4659              } else if (field=="EOL") {
4660                parameters[iParam].printString();
4661                break;
4662              } else {
4663                parameters[iParam].setStringValue(field);
4664              }
4665              std::string fileName;
4666              bool canOpen=false;
4667              if (field=="-") {
4668                // stdin
4669                canOpen=true;
4670                fileName = "-";
4671              } else {
4672                bool absolutePath;
4673                if (dirsep=='/') {
4674                  // non Windows (or cygwin)
4675                  absolutePath=(field[0]=='/');
4676                } else {
4677                  //Windows (non cycgwin)
4678                  absolutePath=(field[0]=='\\');
4679                  // but allow for :
4680                  if (strchr(field.c_str(),':'))
4681                    absolutePath=true;
4682                }
4683                if (absolutePath) {
4684                  fileName = field;
4685                } else if (field[0]=='~') {
4686                  char * environVar = getenv("HOME");
4687                  if (environVar) {
4688                    std::string home(environVar);
4689                    field=field.erase(0,1);
4690                    fileName = home+field;
4691                  } else {
4692                    fileName=field;
4693                  }
4694                } else {
4695                  fileName = directory+field;
4696                }
4697                FILE *fp=fopen(fileName.c_str(),"r");
4698                if (fp) {
4699                  // can open - lets go for it
4700                  fclose(fp);
4701                  canOpen=true;
4702                } else {
4703                  std::cout<<"Unable to open file "<<fileName<<std::endl;
4704                }
4705              }
4706              if (canOpen) {
4707                CoinModel coinModel(fileName.c_str(),2);
4708                // load from coin model
4709                OsiSolverLink solver1;
4710                OsiSolverInterface * solver2 = solver1.clone();
4711                model.assignSolver(solver2,true);
4712                OsiSolverLink * si =
4713                  dynamic_cast<OsiSolverLink *>(model.solver()) ;
4714                assert (si != NULL);
4715                si->setDefaultMeshSize(0.001);
4716                // need some relative granularity
4717                si->setDefaultBound(100.0);
4718                si->setDefaultMeshSize(0.01);
4719                si->setDefaultBound(100.0);
4720                si->setIntegerPriority(1000);
4721                si->setBiLinearPriority(10000);
4722                CoinModel * model2 = (CoinModel *) &coinModel;
4723                si->load(*model2);
4724                // redo
4725                solver = model.solver();
4726                clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
4727                lpSolver = clpSolver->getModelPtr();
4728                clpSolver->messageHandler()->setLogLevel(0) ;
4729                testOsiParameters=0;
4730                complicatedInteger=2;
4731              }
4732            }
4733#endif
4734            break;
4735          case EXPORT:
4736            if (goodModel) {
4737              // get next field
4738              field = CoinReadGetString(argc,argv);
4739              if (field=="$") {
4740                field = parameters[iParam].stringValue();
4741              } else if (field=="EOL") {
4742                parameters[iParam].printString();
4743                break;
4744              } else {
4745                parameters[iParam].setStringValue(field);
4746              }
4747              std::string fileName;
4748              bool canOpen=false;
4749              if (field[0]=='/'||field[0]=='\\') {
4750                fileName = field;
4751              } else if (field[0]=='~') {
4752                char * environVar = getenv("HOME");
4753                if (environVar) {
4754                  std::string home(environVar);
4755                  field=field.erase(0,1);
4756                  fileName = home+field;
4757                } else {
4758                  fileName=field;
4759                }
4760              } else {
4761                fileName = directory+field;
4762              }
4763              FILE *fp=fopen(fileName.c_str(),"w");
4764              if (fp) {
4765                // can open - lets go for it
4766                fclose(fp);
4767                canOpen=true;
4768              } else {
4769                std::cout<<"Unable to open file "<<fileName<<std::endl;
4770              }
4771              if (canOpen) {
4772                // If presolve on then save presolved
4773                bool deleteModel2=false;
4774                ClpSimplex * model2 = lpSolver;
4775                if (dualize) {
4776                  model2 = ((ClpSimplexOther *) model2)->dualOfModel();
4777                  printf("Dual of model has %d rows and %d columns\n",
4778                         model2->numberRows(),model2->numberColumns());
4779                  model2->setOptimizationDirection(1.0);
4780                }
4781#ifdef COIN_HAS_ASL
4782                if (info.numberSos&&doSOS&&usingAmpl) {
4783                  // SOS
4784                  numberSOS = info.numberSos;
4785                  sosStart = info.sosStart;
4786                  sosIndices = info.sosIndices;
4787                  sosReference = info.sosReference;
4788                  preSolve=false;
4789                  clpSolver->setSOSData(numberSOS,info.sosType,sosStart,sosIndices,sosReference);
4790                }
4791#endif
4792                if (preSolve) {
4793                  ClpPresolve pinfo;
4794                  int presolveOptions2 = presolveOptions&~0x40000000;
4795                  if ((presolveOptions2&0xffff)!=0)
4796                    pinfo.setPresolveActions(presolveOptions2);
4797                  if ((printOptions&1)!=0)
4798                    pinfo.statistics();
4799                  double presolveTolerance = 
4800                    parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
4801                  model2 = 
4802                    pinfo.presolvedModel(*lpSolver,presolveTolerance,
4803                                         true,preSolve);
4804                  if (model2) {
4805                    printf("Saving presolved model on %s\n",
4806                           fileName.c_str());
4807                    deleteModel2=true;
4808                  } else {
4809                    printf("Presolved model looks infeasible - saving original on %s\n",
4810                           fileName.c_str());
4811                    deleteModel2=false;
4812                    model2 = lpSolver;
4813
4814                  }
4815                  model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
4816                  if (deleteModel2)
4817                    delete model2;
4818                } else {
4819                  printf("Saving model on %s\n",
4820                           fileName.c_str());
4821                  if (numberSOS) {
4822                    // Convert names
4823                    int iRow;
4824                    int numberRows=model2->numberRows();
4825                    int iColumn;
4826                    int numberColumns=model2->numberColumns();
4827                   
4828                    char ** rowNames = NULL;
4829                    char ** columnNames = NULL;
4830                    if (model2->lengthNames()) {
4831                      rowNames = new char * [numberRows];
4832                      for (iRow=0;iRow<numberRows;iRow++) {
4833                        rowNames[iRow] = 
4834                          strdup(model2->rowName(iRow).c_str());
4835                      }
4836                     
4837                      columnNames = new char * [numberColumns];
4838                      for (iColumn=0;iColumn<numberColumns;iColumn++) {
4839                        columnNames[iColumn] = 
4840                          strdup(model2->columnName(iColumn).c_str());
4841                      }
4842                    }
4843                    clpSolver->writeMpsNative(fileName.c_str(),(const char **) rowNames,(const char **) columnNames,
4844                                              (outputFormat-1)/2,1+((outputFormat-1)&1));
4845                    if (rowNames) {
4846                      for (iRow=0;iRow<numberRows;iRow++) {
4847                        free(rowNames[iRow]);
4848                      }
4849                      delete [] rowNames;
4850                      for (iColumn=0;iColumn<numberColumns;iColumn++) {
4851                        free(columnNames[iColumn]);
4852                      }
4853                      delete [] columnNames;
4854                    }
4855                  } else {
4856                    model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
4857                  }
4858                }
4859                time2 = CoinCpuTime();
4860                totalTime += time2-time1;
4861                time1=time2;
4862              }
4863            } else {
4864              std::cout<<"** Current model not valid"<<std::endl;
4865            }
4866            break;
4867          case BASISIN:
4868            if (goodModel) {
4869              // get next field
4870              field = CoinReadGetString(argc,argv);
4871              if (field=="$") {
4872                field = parameters[iParam].stringValue();
4873              } else if (field=="EOL") {
4874                parameters[iParam].printString();
4875                break;
4876              } else {
4877                parameters[iParam].setStringValue(field);
4878              }
4879              std::string fileName;
4880              bool canOpen=false;
4881              if (field=="-") {
4882                // stdin
4883                canOpen=true;
4884                fileName = "-";
4885              } else {
4886                if (field[0]=='/'||field[0]=='\\') {
4887                  fileName = field;
4888                } else if (field[0]=='~') {
4889                  char * environVar = getenv("HOME");
4890                  if (environVar) {
4891                    std::string home(environVar);
4892                    field=field.erase(0,1);
4893                    fileName = home+field;
4894                  } else {
4895                    fileName=field;
4896                  }
4897                } else {
4898                  fileName = directory+field;
4899                }
4900                FILE *fp=fopen(fileName.c_str(),"r");
4901                if (fp) {
4902                  // can open - lets go for it
4903                  fclose(fp);
4904                  canOpen=true;
4905                } else {
4906                  std::cout<<"Unable to open file "<<fileName<<std::endl;
4907                }
4908              }
4909              if (canOpen) {
4910                int values = lpSolver->readBasis(fileName.c_str());
4911                if (values==0)
4912                  basisHasValues=-1;
4913                else
4914                  basisHasValues=1;
4915              }
4916            } else {
4917              std::cout<<"** Current model not valid"<<std::endl;
4918            }
4919            break;
4920          case PRIORITYIN:
4921            if (goodModel) {
4922              // get next field
4923              field = CoinReadGetString(argc,argv);
4924              if (field=="$") {
4925                field = parameters[iParam].stringValue();
4926              } else if (field=="EOL") {
4927                parameters[iParam].printString();
4928                break;
4929              } else {
4930                parameters[iParam].setStringValue(field);
4931              }
4932              std::string fileName;
4933              if (field[0]=='/'||field[0]=='\\') {
4934                fileName = field;
4935              } else if (field[0]=='~') {
4936                char * environVar = getenv("HOME");
4937                if (environVar) {
4938                  std::string home(environVar);
4939                  field=field.erase(0,1);
4940                  fileName = home+field;
4941                } else {
4942                  fileName=field;
4943                }
4944              } else {
4945                fileName = directory+field;
4946              }
4947              FILE *fp=fopen(fileName.c_str(),"r");
4948              if (fp) {
4949                // can open - lets go for it
4950                std::string headings[]={"name","number","direction","priority","up","down",
4951                                        "solution","priin"};
4952                int got[]={-1,-1,-1,-1,-1,-1,-1,-1};
4953                int order[8];
4954                assert(sizeof(got)==sizeof(order));
4955                int nAcross=0;
4956                char line[1000];
4957                int numberColumns = lpSolver->numberColumns();
4958                if (!fgets(line,1000,fp)) {
4959                  std::cout<<"Odd file "<<fileName<<std::endl;
4960                } else {
4961                  char * pos = line;
4962                  char * put = line;
4963                  while (*pos>=' '&&*pos!='\n') {
4964                    if (*pos!=' '&&*pos!='\t') {
4965                      *put=tolower(*pos);
4966                      put++;
4967                    }
4968                    pos++;
4969                  }
4970                  *put='\0';
4971                  pos=line;
4972                  int i;
4973                  bool good=true;
4974                  while (pos) {
4975                    char * comma = strchr(pos,',');
4976                    if (comma)
4977                      *comma='\0';
4978                    for (i=0;i<(int) (sizeof(got)/sizeof(int));i++) {
4979                      if (headings[i]==pos) {
4980                        if (got[i]<0) {
4981                          order[nAcross]=i;
4982                          got[i]=nAcross++;
4983                        } else {
4984                          // duplicate
4985                          good=false;
4986                        }
4987                        break;
4988                      }
4989                    }
4990                    if (i==(int) (sizeof(got)/sizeof(int)))
4991                      good=false;
4992                    if (comma) {
4993                      *comma=',';
4994                      pos=comma+1;
4995                    } else {
4996                      break;
4997                    }
4998                  }
4999                  if (got[0]<0&&got[1]<0)
5000                    good=false;
5001                  if (got[0]>=0&&got[1]>=0)
5002                    good=false;
5003                  if (got[0]>=0&&!lpSolver->lengthNames())
5004                    good=false;
5005                  if (good) {
5006                    char ** columnNames = new