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

Last change on this file since 679 was 679, checked in by forrest, 14 years ago

take out misleading setting of integer variables because of ampl

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