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

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

out link while I think

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