source: trunk/Cbc/src/CbcSolver.cpp @ 645

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

version

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