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

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

changes for postprocess loop and LOS

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