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

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

for dubiois optionas

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