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

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

message handling

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