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

Last change on this file since 589 was 589, checked in by forrest, 14 years ago

add returncode

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