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

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

nonlinear

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