source: trunk/Cbc/src/CbcSolver.cpp @ 720

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

message handling and sos to CbcSolver?

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