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

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

to make easier to use as function

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