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

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

undef fae

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