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

Last change on this file since 648 was 648, checked in by forrest, 12 years ago

for parallel

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