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

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

changes for gcc 4.3

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