source: branches/BSP/trunk/Cbc/src/CbcSolver.cpp @ 733

Last change on this file since 733 was 733, checked in by andreasw, 14 years ago

corrected Cbc version number that is printed

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