source: branches/devel/Cbc/src/CbcSolver.cpp @ 674

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

timing

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