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

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

for tighter linear constraints in bilinear

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