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

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

ampl mistake

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