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

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

max times and fix bug in fpump

File size: 250.8 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                babModel->setOriginalColumns(process.originalColumns());
3412#ifdef COIN_HAS_ASL
3413                if (!solver2&&usingAmpl) {
3414                  // infeasible
3415                  info.problemStatus=1;
3416                  info.objValue = 1.0e100;
3417                  sprintf(info.buffer,"infeasible/unbounded by pre-processing");
3418                  info.primalSolution=NULL;
3419                  info.dualSolution=NULL;
3420                  break;
3421                }
3422#endif
3423                if (!noPrinting) {
3424                  if (!solver2) {
3425                    sprintf(generalPrint,"Pre-processing says infeasible or unbounded\n");
3426                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
3427                      << generalPrint
3428                      <<CoinMessageEol;
3429                  } else {
3430                    //printf("processed model has %d rows, %d columns and %d elements\n",
3431                    //     solver2->getNumRows(),solver2->getNumCols(),solver2->getNumElements());
3432                  }
3433                }
3434                if (!solver2) {
3435                  model.setProblemStatus(0);
3436                  model.setSecondaryStatus(1);
3437                  babModel->setProblemStatus(0);
3438                  babModel->setSecondaryStatus(1);
3439                } else {
3440                  model.setProblemStatus(-1);
3441                  babModel->setProblemStatus(-1);
3442                }
3443                int returnCode=callBack(babModel,2);
3444                if (returnCode) {
3445                  // exit if user wants
3446                  delete babModel;
3447                  return returnCode;
3448                }
3449                if (!solver2)
3450                  break;
3451                //solver2->resolve();
3452                if (preProcess==2) {
3453                  OsiClpSolverInterface * clpSolver2 = dynamic_cast< OsiClpSolverInterface*> (solver2);
3454                  ClpSimplex * lpSolver = clpSolver2->getModelPtr();
3455                  lpSolver->writeMps("presolved.mps",0,1,lpSolver->optimizationDirection());
3456                  printf("Preprocessed model (minimization) on presolved.mps\n");
3457                }
3458                // we have to keep solver2 so pass clone
3459                solver2 = solver2->clone();
3460                babModel->assignSolver(solver2);
3461                babModel->initialSolve();
3462                babModel->setMaximumSeconds(timeLeft-(CoinCpuTime()-time1));
3463              }
3464              // now tighten bounds
3465              if (!miplib) {
3466                OsiClpSolverInterface * si =
3467                  dynamic_cast<OsiClpSolverInterface *>(babModel->solver()) ;
3468                assert (si != NULL);
3469                // get clp itself
3470                ClpSimplex * modelC = si->getModelPtr();
3471                //if (noPrinting)
3472                //modelC->setLogLevel(0);
3473                if (!complicatedInteger&&modelC->tightenPrimalBounds()!=0) {
3474                  std::cout<<"Problem is infeasible!"<<std::endl;
3475                  break;
3476                }
3477                si->resolve();
3478              }
3479#if 0
3480              numberDebugValues=599;
3481              debugValues = new double[numberDebugValues];
3482              CoinZeroN(debugValues,numberDebugValues);
3483              debugValues[3]=1.0;
3484              debugValues[6]=25.0;
3485              debugValues[9]=4.0;
3486              debugValues[26]=4.0;
3487              debugValues[27]=6.0;
3488              debugValues[35]=8.0;
3489              debugValues[53]=21.0;
3490              debugValues[56]=4.0;
3491#endif
3492              if (debugValues) {
3493                // for debug
3494                std::string problemName ;
3495                babModel->solver()->getStrParam(OsiProbName,problemName) ;
3496                babModel->solver()->activateRowCutDebugger(problemName.c_str()) ;
3497                twomirGen.probname_=strdup(problemName.c_str());
3498                // checking seems odd
3499                //redsplitGen.set_given_optsol(babModel->solver()->getRowCutDebuggerAlways()->optimalSolution(),
3500                //                         babModel->getNumCols());
3501              }
3502              int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
3503#ifdef COIN_HAS_ASL
3504              // If linked then see if expansion wanted
3505              {
3506                OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
3507                int options = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()/10000;
3508                if (solver3||(options&16)!=0) {
3509                  if (options) {
3510                    /*
3511                      1 - force mini branch and bound
3512                      2 - set priorities high on continuous
3513                      4 - try adding OA cuts
3514                      8 - try doing quadratic linearization
3515                      16 - try expanding knapsacks
3516                    */
3517                    if ((options&16)) {
3518                      int numberColumns = saveCoinModel.numberColumns();
3519                      int numberRows = saveCoinModel.numberRows();
3520                      whichColumn = new int[numberColumns];
3521                      knapsackStart=new int[numberRows+1];
3522                      knapsackRow=new int[numberRows];
3523                      numberKnapsack=10000;
3524                      int extra1 = parameters[whichParam(EXTRA1,numberParameters,parameters)].intValue();
3525                      int extra2 = parameters[whichParam(EXTRA2,numberParameters,parameters)].intValue();
3526                      int logLevel = parameters[log].intValue();
3527                      OsiSolverInterface * solver = expandKnapsack(saveCoinModel,whichColumn,knapsackStart,
3528                                                                   knapsackRow,numberKnapsack,
3529                                                                   storedAmpl,logLevel,extra1,extra2,
3530                                                                   saveTightenedModel);
3531                      if (solver) {
3532                        clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3533                        assert (clpSolver);
3534                        lpSolver = clpSolver->getModelPtr();
3535                        babModel->assignSolver(solver);
3536                        testOsiOptions=0;
3537                        // allow gomory
3538                        complicatedInteger=0;
3539                        // Priorities already done
3540                        free(info.priorities);
3541                        info.priorities=NULL;
3542                      } else {
3543                        numberKnapsack=0;
3544                        delete [] whichColumn;
3545                        delete [] knapsackStart;
3546                        delete [] knapsackRow;
3547                        whichColumn=NULL;
3548                        knapsackStart=NULL;
3549                        knapsackRow=NULL;
3550                      }
3551                    }
3552                  }
3553                }
3554              }
3555#endif
3556              if (useCosts&&testOsiOptions<0) {
3557                int numberColumns = babModel->getNumCols();
3558                int * sort = new int[numberColumns];
3559                double * dsort = new double[numberColumns];
3560                int * priority = new int [numberColumns];
3561                const double * objective = babModel->getObjCoefficients();
3562                const double * lower = babModel->getColLower() ;
3563                const double * upper = babModel->getColUpper() ;
3564                const CoinPackedMatrix * matrix = babModel->solver()->getMatrixByCol();
3565                const int * columnLength = matrix->getVectorLengths();
3566                int iColumn;
3567                int n=0;
3568                for (iColumn=0;iColumn<numberColumns;iColumn++) {
3569                  if (babModel->isInteger(iColumn)) {
3570                    sort[n]=n;
3571                    if (useCosts==1)
3572                      dsort[n++]=-fabs(objective[iColumn]);
3573                    else if (useCosts==2)
3574                      dsort[n++]=iColumn;
3575                    else if (useCosts==3)
3576                      dsort[n++]=upper[iColumn]-lower[iColumn];
3577                    else if (useCosts==4)
3578                      dsort[n++]=-(upper[iColumn]-lower[iColumn]);
3579                    else if (useCosts==5)
3580                      dsort[n++]=-columnLength[iColumn];
3581                  }
3582                }
3583                CoinSort_2(dsort,dsort+n,sort);
3584                int level=0;
3585                double last = -1.0e100;
3586                for (int i=0;i<n;i++) {
3587                  int iPut=sort[i];
3588                  if (dsort[i]!=last) {
3589                    level++;
3590                    last=dsort[i];
3591                  }
3592                  priority[iPut]=level;
3593                }
3594                babModel->passInPriorities( priority,false);
3595                integersOK=true;
3596                delete [] priority;
3597                delete [] sort;
3598                delete [] dsort;
3599              }
3600              // FPump done first as it only works if no solution
3601              CbcHeuristicFPump heuristic4(*babModel);
3602              heuristic4.setFractionSmall(0.6);
3603              if (useFpump) {
3604                heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
3605                int pumpTune=parameters[whichParam(FPUMPTUNE,numberParameters,parameters)].intValue();
3606                if (pumpTune>0) {
3607                  /*
3608                    >=10000000 for using obj
3609                    >=1000000 use as accumulate switch
3610                    >=1000 use index+1 as number of large loops
3611                    >=100 use 0.05 objvalue as increment
3612                    >=10 use +0.1 objvalue for cutoff (add)
3613                    1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
3614                    4 and static continuous, 5 as 3 but no internal integers
3615                    6 as 3 but all slack basis!
3616                  */
3617                  int logLevel = parameters[log].intValue();
3618                  double value = babModel->solver()->getObjSense()*babModel->solver()->getObjValue();
3619                  int w = pumpTune/10;
3620                  int c = w % 10;
3621                  w /= 10;
3622                  int i = w % 10;
3623                  w /= 10;
3624                  int r = w;
3625                  int accumulate = r/1000;
3626                  r -= 1000*accumulate;
3627                  if (accumulate>=10) {
3628                    int which = accumulate/10;
3629                    accumulate -= 10*which;
3630                    which--;
3631                    // weights and factors
3632                    double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
3633                    double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
3634                    heuristic4.setInitialWeight(weight[which]);
3635                    heuristic4.setWeightFactor(factor[which]);
3636                  }
3637                  // fake cutoff
3638                  if (logLevel>1)
3639                    printf("Setting ");
3640                  if (c) {
3641                    double cutoff;
3642                    babModel->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
3643                    cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
3644                    double dextra1 = parameters[whichParam(DEXTRA1,numberParameters,parameters)].doubleValue();
3645                    if (dextra1)
3646                      cutoff=dextra1;
3647                    heuristic4.setFakeCutoff(cutoff);
3648                    if (logLevel>1)
3649                      printf("fake cutoff of %g ",cutoff);
3650                  }
3651                  if (i||r) {
3652                    // also set increment
3653                    double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
3654                    double dextra2 = parameters[whichParam(DEXTRA2,numberParameters,parameters)].doubleValue();
3655                    if (dextra2)
3656                      increment = dextra2;
3657                    heuristic4.setAbsoluteIncrement(increment);
3658                    heuristic4.setAccumulate(accumulate);
3659                    heuristic4.setMaximumRetries(r+1);
3660                    if (logLevel>1) {
3661                      if (i) 
3662                        printf("increment of %g ",heuristic4.absoluteIncrement());
3663                      if (accumulate) 
3664                        printf("accumulate of %d ",accumulate);
3665                      printf("%d retries ",r+2);
3666                    }
3667                  }
3668                  pumpTune = pumpTune%100;
3669                  if (logLevel>1)
3670                    printf("and setting when to %d\n",pumpTune+10);
3671                  if (pumpTune==6)
3672                    pumpTune =13;
3673                  heuristic4.setWhen(pumpTune+10);
3674                }
3675                heuristic4.setHeuristicName("feasibility pump");
3676                babModel->addHeuristic(&heuristic4);
3677              }
3678              if (!miplib) {
3679                CbcRounding heuristic1(*babModel);
3680                heuristic1.setHeuristicName("rounding");
3681                if (useRounding)
3682                  babModel->addHeuristic(&heuristic1) ;
3683                CbcHeuristicLocal heuristic2(*babModel);
3684                heuristic2.setHeuristicName("combine solutions");
3685                heuristic2.setFractionSmall(0.6);
3686                heuristic2.setSearchType(1);
3687                if (useCombine)
3688                  babModel->addHeuristic(&heuristic2);
3689                CbcHeuristicGreedyCover heuristic3(*babModel);
3690                heuristic3.setHeuristicName("greedy cover");
3691                CbcHeuristicGreedyEquality heuristic3a(*babModel);
3692                heuristic3a.setHeuristicName("greedy equality");
3693                if (useGreedy) {
3694                  babModel->addHeuristic(&heuristic3);
3695                  babModel->addHeuristic(&heuristic3a);
3696                }
3697                if (useLocalTree) {
3698                  CbcTreeLocal localTree(babModel,NULL,10,0,0,10000,2000);
3699                  babModel->passInTreeHandler(localTree);
3700                }
3701              }
3702              CbcHeuristicRINS heuristic5(*babModel);
3703              heuristic5.setHeuristicName("RINS");
3704              heuristic5.setFractionSmall(0.6);
3705              if (useRINS)
3706                babModel->addHeuristic(&heuristic5) ;
3707              if (type==MIPLIB) {
3708                if (babModel->numberStrong()==5&&babModel->numberBeforeTrust()==5) 
3709                  babModel->setNumberBeforeTrust(10);
3710              }
3711              // add cut generators if wanted
3712              int switches[20];
3713              int numberGenerators=0;
3714              int translate[]={-100,-1,-99,-98,1,1,1,1};
3715              if (probingAction) {
3716                if (probingAction==5||probingAction==7)
3717                  probingGen.setRowCuts(-3); // strengthening etc just at root
3718                if (probingAction==6||probingAction==7) {
3719                  // Number of unsatisfied variables to look at
3720                  probingGen.setMaxProbe(1000);
3721                  probingGen.setMaxProbeRoot(1000);
3722                  // How far to follow the consequences
3723                  probingGen.setMaxLook(50);
3724                  probingGen.setMaxLookRoot(50);
3725                }
3726                babModel->addCutGenerator(&probingGen,translate[probingAction],"Probing");
3727                switches[numberGenerators++]=0;
3728              }
3729              if (gomoryAction&&(complicatedInteger!=1||
3730                                 (gomoryAction==1||gomoryAction==4))) {
3731                babModel->addCutGenerator(&gomoryGen,translate[gomoryAction],"Gomory");
3732                switches[numberGenerators++]=-1;
3733              }
3734              if (knapsackAction) {
3735                babModel->addCutGenerator(&knapsackGen,translate[knapsackAction],"Knapsack");
3736                switches[numberGenerators++]=0;
3737              }
3738              if (redsplitAction&&!complicatedInteger) {
3739                babModel->addCutGenerator(&redsplitGen,translate[redsplitAction],"Reduce-and-split");
3740                switches[numberGenerators++]=1;
3741              }
3742              if (cliqueAction) {
3743                babModel->addCutGenerator(&cliqueGen,translate[cliqueAction],"Clique");
3744                switches[numberGenerators++]=0;
3745              }
3746              if (mixedAction) {
3747                babModel->addCutGenerator(&mixedGen,translate[mixedAction],"MixedIntegerRounding2");
3748                switches[numberGenerators++]=-1;
3749              }
3750              if (flowAction) {
3751                babModel->addCutGenerator(&flowGen,translate[flowAction],"FlowCover");
3752                switches[numberGenerators++]=1;
3753              }
3754              if (twomirAction&&!complicatedInteger) {
3755                babModel->addCutGenerator(&twomirGen,translate[twomirAction],"TwoMirCuts");
3756                switches[numberGenerators++]=1;
3757              }
3758              if (landpAction) {
3759                babModel->addCutGenerator(&landpGen,translate[landpAction],"LiftAndProject");
3760                switches[numberGenerators++]=1;
3761              }
3762              if (residualCapacityAction) {
3763                babModel->addCutGenerator(&residualCapacityGen,translate[residualCapacityAction],"ResidualCapacity");
3764                switches[numberGenerators++]=1;
3765              }
3766              if (storedCuts) 
3767                babModel->setSpecialOptions(babModel->specialOptions()|64);
3768              // Say we want timings
3769              numberGenerators = babModel->numberCutGenerators();
3770              int iGenerator;
3771              int cutDepth=
3772                parameters[whichParam(CUTDEPTH,numberParameters,parameters)].intValue();
3773              for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
3774                CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
3775                int howOften = generator->howOften();
3776                if (howOften==-98||howOften==-99) 
3777                  generator->setSwitchOffIfLessThan(switches[iGenerator]);
3778                generator->setTiming(true);
3779                if (cutDepth>=0)
3780                  generator->setWhatDepth(cutDepth) ;
3781              }
3782              // Could tune more
3783              if (!miplib) {
3784                babModel->setMinimumDrop(min(5.0e-2,
3785                                             fabs(babModel->getMinimizationObjValue())*1.0e-3+1.0e-4));
3786                if (cutPass==-1234567) {
3787                  if (babModel->getNumCols()<500)
3788                    babModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
3789                  else if (babModel->getNumCols()<5000)
3790                    babModel->setMaximumCutPassesAtRoot(100); // use minimum drop
3791                  else
3792                    babModel->setMaximumCutPassesAtRoot(20);
3793                } else {
3794                  babModel->setMaximumCutPassesAtRoot(cutPass);
3795                }
3796                if (cutPassInTree==-1234567) 
3797                  babModel->setMaximumCutPasses(1);
3798                else
3799                  babModel->setMaximumCutPasses(cutPassInTree);
3800              }
3801              // Do more strong branching if small
3802              //if (babModel->getNumCols()<5000)
3803              //babModel->setNumberStrong(20);
3804              // Switch off strong branching if wanted
3805              //if (babModel->getNumCols()>10*babModel->getNumRows())
3806              //babModel->setNumberStrong(0);
3807              if (!noPrinting) {
3808                int iLevel = parameters[log].intValue();
3809                if (iLevel<0) {
3810                  babModel->setPrintingMode(1);
3811                  iLevel = -iLevel;
3812                }
3813                babModel->messageHandler()->setLogLevel(iLevel);
3814                if (babModel->getNumCols()>2000||babModel->getNumRows()>1500||
3815                    babModel->messageHandler()->logLevel()>1)
3816                  babModel->setPrintFrequency(100);
3817              }
3818             
3819              babModel->solver()->setIntParam(OsiMaxNumIterationHotStart,
3820                    parameters[whichParam(MAXHOTITS,numberParameters,parameters)].intValue());
3821              OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
3822              // go faster stripes
3823              if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
3824                osiclp->setupForRepeatedUse(2,parameters[slog].intValue());
3825              } else {
3826                osiclp->setupForRepeatedUse(0,parameters[slog].intValue());
3827              }
3828              double increment=babModel->getCutoffIncrement();;
3829              int * changed = NULL;
3830              if (!miplib&&increment==normalIncrement)
3831                changed=analyze( osiclp,numberChanged,increment,false,generalMessageHandler);
3832              if (debugValues) {
3833                if (numberDebugValues==babModel->getNumCols()) {
3834                  // for debug
3835                  babModel->solver()->activateRowCutDebugger(debugValues) ;
3836                } else {
3837                  printf("debug file has incorrect number of columns\n");
3838                }
3839              }
3840              babModel->setCutoffIncrement(CoinMax(babModel->getCutoffIncrement(),increment));
3841              // Turn this off if you get problems
3842              // Used to be automatically set
3843              int mipOptions = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()%10000;
3844              if (mipOptions!=(1)) {
3845                sprintf(generalPrint,"mip options %d\n",mipOptions);
3846                generalMessageHandler->message(CLP_GENERAL,generalMessages)
3847                  << generalPrint
3848                  <<CoinMessageEol;
3849              }
3850              osiclp->setSpecialOptions(mipOptions);
3851              if (gapRatio < 1.0e100) {
3852                double value = babModel->solver()->getObjValue() ;
3853                double value2 = gapRatio*(1.0e-5+fabs(value)) ;
3854                babModel->setAllowableGap(value2) ;
3855                std::cout << "Continuous " << value
3856                          << ", so allowable gap set to "
3857                          << value2 << std::endl ;
3858              }
3859              // probably faster to use a basis to get integer solutions
3860              babModel->setSpecialOptions(babModel->specialOptions()|2);
3861              currentBranchModel = babModel;
3862              OsiSolverInterface * strengthenedModel=NULL;
3863              if (type==BAB||type==MIPLIB) {
3864                int moreMipOptions = parameters[whichParam(MOREMIPOPTIONS,numberParameters,parameters)].intValue();
3865                if (moreMipOptions>=0) {
3866                  sprintf(generalPrint,"more mip options %d\n",moreMipOptions);
3867                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
3868                    << generalPrint
3869                    <<CoinMessageEol;
3870                  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
3871                  if (moreMipOptions==10000) {
3872                    // test memory saving
3873                    moreMipOptions -= 10000;
3874                    ClpSimplex * lpSolver = osiclp->getModelPtr();
3875                    lpSolver->setPersistenceFlag(1);
3876                    // switch off row copy if few rows
3877                    if (lpSolver->numberRows()<150)
3878                      lpSolver->setSpecialOptions(lpSolver->specialOptions()|256);
3879                  }
3880                  if (((moreMipOptions+1)%1000000)!=0)
3881                    babModel->setSearchStrategy(moreMipOptions%1000000);
3882                  // go faster stripes
3883                  if( moreMipOptions >=999999) {
3884                    if (osiclp) {
3885                      int save = osiclp->specialOptions();
3886                      osiclp->setupForRepeatedUse(2,0);
3887                      osiclp->setSpecialOptions(save|osiclp->specialOptions());
3888                    }
3889                  } 
3890                }
3891              }
3892              if (type==BAB) {
3893#ifdef COIN_HAS_ASL
3894                if (usingAmpl) {
3895                  priorities=info.priorities;
3896                  branchDirection=info.branchDirection;
3897                  pseudoDown=info.pseudoDown;
3898                  pseudoUp=info.pseudoUp;
3899                  solutionIn=info.primalSolution;
3900                  prioritiesIn = info.priorities;
3901                  if (info.numberSos&&doSOS) {
3902                    // SOS
3903                    numberSOS = info.numberSos;
3904                    sosStart = info.sosStart;
3905                    sosIndices = info.sosIndices;
3906                    sosType = info.sosType;
3907                    sosReference = info.sosReference;
3908                    sosPriority = info.sosPriority;
3909                  }
3910                }
3911#endif               
3912                const int * originalColumns = preProcess ? process.originalColumns() : NULL;
3913                if (solutionIn&&useSolution) {
3914                  if (preProcess) {
3915                    int numberColumns = babModel->getNumCols();
3916                    // extend arrays in case SOS
3917                    int n = originalColumns[numberColumns-1]+1;
3918                    int nSmaller = CoinMin(n,numberOriginalColumns);
3919                    double * solutionIn2 = new double [n];
3920                    int * prioritiesIn2 = new int[n];
3921                    int i;
3922                    for (i=0;i<nSmaller;i++) {
3923                      solutionIn2[i]=solutionIn[i];
3924                      prioritiesIn2[i]=prioritiesIn[i];
3925                    }
3926                    for (;i<n;i++) {
3927                      solutionIn2[i]=0.0;
3928                      prioritiesIn2[i]=1000000;
3929                    }
3930                    int iLast=-1;
3931                    for (i=0;i<numberColumns;i++) {
3932                      int iColumn = originalColumns[i];
3933                      assert (iColumn>iLast);
3934                      iLast=iColumn;
3935                      solutionIn2[i]=solutionIn2[iColumn];
3936                      if (prioritiesIn)
3937                        prioritiesIn2[i]=prioritiesIn2[iColumn];
3938                    }
3939                    babModel->setHotstartSolution(solutionIn2,prioritiesIn2);
3940                    delete [] solutionIn2;
3941                    delete [] prioritiesIn2;
3942                  } else {
3943                    babModel->setHotstartSolution(solutionIn,prioritiesIn);
3944                  }
3945                }
3946                OsiSolverInterface * testOsiSolver= (testOsiOptions>=0) ? babModel->solver() : NULL;
3947                if (!testOsiSolver) {
3948                  // *************************************************************
3949                  // CbcObjects
3950                  if (preProcess&&(process.numberSOS()||babModel->numberObjects())) {
3951                    int numberSOS = process.numberSOS();
3952                    int numberIntegers = babModel->numberIntegers();
3953                    /* model may not have created objects
3954                       If none then create
3955                    */
3956                    if (!numberIntegers||!babModel->numberObjects()) {
3957                      int type = (pseudoUp) ? 1 : 0;
3958                      babModel->findIntegers(true,type);
3959                      numberIntegers = babModel->numberIntegers();
3960                      integersOK=true;
3961                    }
3962                    OsiObject ** oldObjects = babModel->objects();
3963                    // Do sets and priorities
3964                    OsiObject ** objects = new OsiObject * [numberSOS];
3965                    // set old objects to have low priority
3966                    int numberOldObjects = babModel->numberObjects();
3967                    int numberColumns = babModel->getNumCols();
3968                    // backward pointer to new variables
3969                    int * newColumn = new int[numberOriginalColumns];
3970                    int i;
3971                    for (i=0;i<numberOriginalColumns;i++)
3972                      newColumn[i]=-1;
3973                    assert (originalColumns);
3974                    for (i=0;i<numberColumns;i++)
3975                      newColumn[originalColumns[i]]=i;
3976                    if (!integersOK) {
3977                      // Change column numbers etc
3978                      int n=0;
3979                      for (int iObj = 0;iObj<numberOldObjects;iObj++) {
3980                        int iColumn = oldObjects[iObj]->columnNumber();
3981                        if (iColumn<0||iColumn>=numberOriginalColumns) {
3982                          oldObjects[n++]=oldObjects[iObj];
3983                        } else {
3984                          iColumn = newColumn[iColumn];
3985                          if (iColumn>=0) {
3986                            CbcSimpleInteger * obj =
3987                              dynamic_cast <CbcSimpleInteger *>(oldObjects[iObj]) ;
3988                            assert (obj);
3989                            obj->setColumnNumber(iColumn);
3990                            oldObjects[n++]=oldObjects[iObj];
3991                          } else {
3992                            delete oldObjects[iObj];
3993                          }
3994                        }
3995                      }
3996                      babModel->setNumberObjects(n);
3997                    }
3998                    int nMissing=0;
3999                    for (int iObj = 0;iObj<numberOldObjects;iObj++) {
4000                      if (process.numberSOS())
4001                        oldObjects[iObj]->setPriority(numberColumns+1);
4002                      int iColumn = oldObjects[iObj]->columnNumber();
4003                      if (iColumn<0||iColumn>=numberOriginalColumns) {
4004                        CbcSOS * obj =
4005                          dynamic_cast <CbcSOS *>(oldObjects[iObj]) ;
4006                        if (obj) {
4007                          int n=obj->numberMembers();
4008                          int * which = obj->mutableMembers();
4009                          double * weights = obj->mutableWeights();
4010                          int nn=0;
4011                          for (i=0;i<n;i++) {
4012                            int iColumn = which[i];
4013                            int jColumn = newColumn[iColumn];
4014                            if (jColumn>=0) { 
4015                              which[nn] = jColumn;
4016                              weights[nn++]=weights[i];
4017                            } else {
4018                              nMissing++;
4019                            }
4020                          }
4021                          obj->setNumberMembers(nn);
4022                        }
4023                        continue;
4024                      }
4025                      if (originalColumns)
4026                        iColumn = originalColumns[iColumn];
4027                      if (branchDirection) {
4028                        CbcSimpleInteger * obj =
4029                          dynamic_cast <CbcSimpleInteger *>(oldObjects[iObj]) ;
4030                        if (obj) { 
4031                          obj->setPreferredWay(branchDirection[iColumn]);
4032                        } else {
4033                          CbcObject * obj =
4034                            dynamic_cast <CbcObject *>(oldObjects[iObj]) ;
4035                          assert (obj);
4036                          obj->setPreferredWay(branchDirection[iColumn]);
4037                        }
4038                      }
4039                      if (pseudoUp) {
4040                        CbcSimpleIntegerPseudoCost * obj1a =
4041                          dynamic_cast <CbcSimpleIntegerPseudoCost *>(oldObjects[iObj]) ;
4042                        assert (obj1a);
4043                        if (pseudoDown[iColumn]>0.0)
4044                          obj1a->setDownPseudoCost(pseudoDown[iColumn]);
4045                        if (pseudoUp[iColumn]>0.0)
4046                          obj1a->setUpPseudoCost(pseudoUp[iColumn]);
4047                      }
4048                    }
4049                    if (nMissing) {
4050                      sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
4051                      generalMessageHandler->message(CLP_GENERAL,generalMessages)
4052                        << generalPrint
4053                        <<CoinMessageEol;
4054                    }
4055                    delete [] newColumn;
4056                    const int * starts = process.startSOS();
4057                    const int * which = process.whichSOS();
4058                    const int * type = process.typeSOS();
4059                    const double * weight = process.weightSOS();
4060                    int iSOS;
4061                    for (iSOS =0;iSOS<numberSOS;iSOS++) {
4062                      int iStart = starts[iSOS];
4063                      int n=starts[iSOS+1]-iStart;
4064                      objects[iSOS] = new CbcSOS(babModel,n,which+iStart,weight+iStart,
4065                                                 iSOS,type[iSOS]);
4066                      // branch on long sets first
4067                      objects[iSOS]->setPriority(numberColumns-n);
4068                    }
4069                    if (numberSOS)
4070                      babModel->addObjects(numberSOS,objects);
4071                    for (iSOS=0;iSOS<numberSOS;iSOS++)
4072                      delete objects[iSOS];
4073                    delete [] objects;
4074                  } else if (priorities||branchDirection||pseudoDown||pseudoUp||numberSOS) {
4075                    // do anyway for priorities etc
4076                    int numberIntegers = babModel->numberIntegers();
4077                    /* model may not have created objects
4078                       If none then create
4079                    */
4080                    if (!numberIntegers||!babModel->numberObjects()) {
4081                      int type = (pseudoUp) ? 1 : 0;
4082                      babModel->findIntegers(true,type);
4083                    }
4084                    if (numberSOS) {
4085                      // Do sets and priorities
4086                      OsiObject ** objects = new OsiObject * [numberSOS];
4087                      int iSOS;
4088                      if (originalColumns) {
4089                        // redo sequence numbers
4090                        int numberColumns = babModel->getNumCols();
4091                        int nOld = originalColumns[numberColumns-1]+1;
4092                        int * back = new int[nOld];
4093                        int i;
4094                        for (i=0;i<nOld;i++)
4095                          back[i]=-1;
4096                        for (i=0;i<numberColumns;i++)
4097                          back[originalColumns[i]]=i;
4098                        // Really need better checks
4099                        int nMissing=0;
4100                        int n=sosStart[numberSOS];
4101                        for (i=0;i<n;i++) {
4102                          int iColumn = sosIndices[i];
4103                          int jColumn = back[iColumn];
4104                          if (jColumn>=0) 
4105                            sosIndices[i] = jColumn;
4106                          else 
4107                            nMissing++;
4108                        }
4109                        delete [] back;
4110                        if (nMissing) {
4111                          sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
4112                          generalMessageHandler->message(CLP_GENERAL,generalMessages)
4113                            << generalPrint
4114                            <<CoinMessageEol;
4115                        }
4116                      }
4117                      for (iSOS =0;iSOS<numberSOS;iSOS++) {
4118                        int iStart = sosStart[iSOS];
4119                        int n=sosStart[iSOS+1]-iStart;
4120                        objects[iSOS] = new CbcSOS(babModel,n,sosIndices+iStart,sosReference+iStart,
4121                                                   iSOS,sosType[iSOS]);
4122                        if (sosPriority)
4123                          objects[iSOS]->setPriority(sosPriority[iSOS]);
4124                        else if (!prioritiesIn)
4125                          objects[iSOS]->setPriority(10);  // rather than 1000
4126                      }
4127                      // delete any existing SOS objects
4128                      int numberObjects=babModel->numberObjects();
4129                      OsiObject ** oldObjects=babModel->objects();
4130                      int nNew=0;
4131                      for (int i=0;i<numberObjects;i++) {
4132                        OsiObject * objThis = oldObjects[i];
4133                        CbcSOS * obj1 =
4134                          dynamic_cast <CbcSOS *>(objThis) ;
4135                        OsiSOS * obj2 =
4136                          dynamic_cast <OsiSOS *>(objThis) ;
4137                        if (!obj1&&!obj2) {
4138                          oldObjects[nNew++]=objThis;
4139                        } else {
4140                          delete objThis;
4141                        }
4142                      }
4143                      babModel->setNumberObjects(nNew);
4144                      babModel->addObjects(numberSOS,objects);
4145                      for (iSOS=0;iSOS<numberSOS;iSOS++)
4146                        delete objects[iSOS];
4147                      delete [] objects;
4148                    }
4149                  }
4150                  OsiObject ** objects = babModel->objects();
4151                  int numberObjects = babModel->numberObjects();
4152                  for (int iObj = 0;iObj<numberObjects;iObj++) {
4153                    // skip sos
4154                    CbcSOS * objSOS =
4155                      dynamic_cast <CbcSOS *>(objects[iObj]) ;
4156                    if (objSOS)
4157                      continue;
4158                    int iColumn = objects[iObj]->columnNumber();
4159                    assert (iColumn>=0);
4160                    if (originalColumns)
4161                      iColumn = originalColumns[iColumn];
4162                    if (branchDirection) {
4163                      CbcSimpleInteger * obj =
4164                        dynamic_cast <CbcSimpleInteger *>(objects[iObj]) ;
4165                      if (obj) { 
4166                        obj->setPreferredWay(branchDirection[iColumn]);
4167                      } else {
4168                        CbcObject * obj =
4169                          dynamic_cast <CbcObject *>(objects[iObj]) ;
4170                        assert (obj);
4171                        obj->setPreferredWay(branchDirection[iColumn]);
4172                      }
4173                    }
4174                    if (priorities) {
4175                      int iPriority = priorities[iColumn];
4176                      if (iPriority>0)
4177                        objects[iObj]->setPriority(iPriority);
4178                    }
4179                    if (pseudoUp&&pseudoUp[iColumn]) {
4180                      CbcSimpleIntegerPseudoCost * obj1a =
4181                        dynamic_cast <CbcSimpleIntegerPseudoCost *>(objects[iObj]) ;
4182                      assert (obj1a);
4183                      if (pseudoDown[iColumn]>0.0)
4184                        obj1a->setDownPseudoCost(pseudoDown[iColumn]);
4185                      if (pseudoUp[iColumn]>0.0)
4186                        obj1a->setUpPseudoCost(pseudoUp[iColumn]);
4187                    }
4188                  }
4189                  // *************************************************************
4190                } else {
4191                  // *************************************************************
4192                  // OsiObjects
4193                  // Find if none
4194                  int numberIntegers = testOsiSolver->getNumIntegers();
4195                  /* model may not have created objects
4196                     If none then create
4197                  */
4198                  if (!numberIntegers||!testOsiSolver->numberObjects()) {
4199                    //int type = (pseudoUp) ? 1 : 0;
4200                    testOsiSolver->findIntegers(false);
4201                    numberIntegers = testOsiSolver->getNumIntegers();
4202                  }
4203                  if (preProcess&&process.numberSOS()) {
4204                    int numberSOS = process.numberSOS();
4205                    OsiObject ** oldObjects = testOsiSolver->objects();
4206                    // Do sets and priorities
4207                    OsiObject ** objects = new OsiObject * [numberSOS];
4208                    // set old objects to have low priority
4209                    int numberOldObjects = testOsiSolver->numberObjects();
4210                    int numberColumns = testOsiSolver->getNumCols();
4211                    for (int iObj = 0;iObj<numberOldObjects;iObj++) {
4212                      oldObjects[iObj]->setPriority(numberColumns+1);
4213                      int iColumn = oldObjects[iObj]->columnNumber();
4214                      assert (iColumn>=0);
4215                      if (iColumn>=numberOriginalColumns)
4216                        continue;
4217                      if (originalColumns)
4218                        iColumn = originalColumns[iColumn];
4219                      if (branchDirection) {
4220                        OsiSimpleInteger * obj =
4221                          dynamic_cast <OsiSimpleInteger *>(oldObjects[iObj]) ;
4222                        if (obj) { 
4223                          obj->setPreferredWay(branchDirection[iColumn]);
4224                        } else {
4225                          OsiObject2 * obj =
4226                            dynamic_cast <OsiObject2 *>(oldObjects[iObj]) ;
4227                          if (obj)
4228                            obj->setPreferredWay(branchDirection[iColumn]);
4229                        }
4230                      }
4231                      if (pseudoUp) {
4232                        abort();
4233                      }
4234                    }
4235                    const int * starts = process.startSOS();
4236                    const int * which = process.whichSOS();
4237                    const int * type = process.typeSOS();
4238                    const double * weight = process.weightSOS();
4239                    int iSOS;
4240                    for (iSOS =0;iSOS<numberSOS;iSOS++) {
4241                      int iStart = starts[iSOS];
4242                      int n=starts[iSOS+1]-iStart;
4243                      objects[iSOS] = new OsiSOS(testOsiSolver,n,which+iStart,weight+iStart,
4244                                                 type[iSOS]);
4245                      // branch on long sets first
4246                      objects[iSOS]->setPriority(numberColumns-n);
4247                    }
4248                    testOsiSolver->addObjects(numberSOS,objects);
4249                    for (iSOS=0;iSOS<numberSOS;iSOS++)
4250                      delete objects[iSOS];
4251                    delete [] objects;
4252                  } else if (priorities||branchDirection||pseudoDown||pseudoUp||numberSOS) {
4253                    if (numberSOS) {
4254                      // Do sets and priorities
4255                      OsiObject ** objects = new OsiObject * [numberSOS];
4256                      int iSOS;
4257                      if (originalColumns) {
4258                        // redo sequence numbers
4259                        int numberColumns = testOsiSolver->getNumCols();
4260                        int nOld = originalColumns[numberColumns-1]+1;
4261                        int * back = new int[nOld];
4262                        int i;
4263                        for (i=0;i<nOld;i++)
4264                          back[i]=-1;
4265                        for (i=0;i<numberColumns;i++)
4266                          back[originalColumns[i]]=i;
4267                        // Really need better checks
4268                        int nMissing=0;
4269                        int n=sosStart[numberSOS];
4270                        for (i=0;i<n;i++) {
4271                          int iColumn = sosIndices[i];
4272                          int jColumn = back[iColumn];
4273                          if (jColumn>=0) 
4274                            sosIndices[i] = jColumn;
4275                          else 
4276                            nMissing++;
4277                        }
4278                        delete [] back;
4279                        if (nMissing) {
4280                          sprintf(generalPrint,"%d SOS variables vanished due to pre processing? - check validity?\n",nMissing);
4281                          generalMessageHandler->message(CLP_GENERAL,generalMessages)
4282                            << generalPrint
4283                            <<CoinMessageEol;
4284                        }
4285                      }
4286                      for (iSOS =0;iSOS<numberSOS;iSOS++) {
4287                        int iStart = sosStart[iSOS];
4288                        int n=sosStart[iSOS+1]-iStart;
4289                        objects[iSOS] = new OsiSOS(testOsiSolver,n,sosIndices+iStart,sosReference+iStart,
4290                                                   sosType[iSOS]);
4291                        if (sosPriority)
4292                          objects[iSOS]->setPriority(sosPriority[iSOS]);
4293                        else if (!prioritiesIn)
4294                          objects[iSOS]->setPriority(10);  // rather than 1000
4295                      }
4296                      // delete any existing SOS objects
4297                      int numberObjects=testOsiSolver->numberObjects();
4298                      OsiObject ** oldObjects=testOsiSolver->objects();
4299                      int nNew=0;
4300                      for (int i=0;i<numberObjects;i++) {
4301                        OsiObject * objThis = oldObjects[i];
4302                        OsiSOS * obj1 =
4303                          dynamic_cast <OsiSOS *>(objThis) ;
4304                        OsiSOS * obj2 =
4305                          dynamic_cast <OsiSOS *>(objThis) ;
4306                        if (!obj1&&!obj2) {
4307                          oldObjects[nNew++]=objThis;
4308                        } else {
4309                          delete objThis;
4310                        }
4311                      }
4312                      testOsiSolver->setNumberObjects(nNew);
4313                      testOsiSolver->addObjects(numberSOS,objects);
4314                      for (iSOS=0;iSOS<numberSOS;iSOS++)
4315                        delete objects[iSOS];
4316                      delete [] objects;
4317                    }
4318                  }
4319                  OsiObject ** objects = testOsiSolver->objects();
4320                  int numberObjects = testOsiSolver->numberObjects();
4321                  int logLevel = parameters[log].intValue();
4322                  for (int iObj = 0;iObj<numberObjects;iObj++) {
4323                    // skip sos
4324                    OsiSOS * objSOS =
4325                      dynamic_cast <OsiSOS *>(objects[iObj]) ;
4326                    if (objSOS) {
4327                      if (logLevel>2)
4328                        printf("Set %d is SOS - priority %d\n",iObj,objSOS->priority());
4329                      continue;
4330                    }
4331                    int iColumn = objects[iObj]->columnNumber();
4332                    if (iColumn>=0) {
4333                      if (originalColumns)
4334                        iColumn = originalColumns[iColumn];
4335                      if (branchDirection) {
4336                        OsiSimpleInteger * obj =
4337                          dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
4338                        if (obj) { 
4339                          obj->setPreferredWay(branchDirection[iColumn]);
4340                        } else {
4341                          OsiObject2 * obj =
4342                            dynamic_cast <OsiObject2 *>(objects[iObj]) ;
4343                          if (obj)
4344                            obj->setPreferredWay(branchDirection[iColumn]);
4345                        }
4346                      }
4347                      if (priorities) {
4348                        int iPriority = priorities[iColumn];
4349                        if (iPriority>0)
4350                          objects[iObj]->setPriority(iPriority);
4351                      }
4352                      if (logLevel>2)
4353                        printf("Obj %d is int? - priority %d\n",iObj,objects[iObj]->priority());
4354                      if (pseudoUp&&pseudoUp[iColumn]) {
4355                        abort();
4356                      }
4357                    }
4358                  }
4359                  // *************************************************************
4360                }
4361                int statistics = (printOptions>0) ? printOptions: 0;
4362#ifdef COIN_HAS_ASL
4363                if (!usingAmpl) {
4364#endif
4365                  free(priorities);
4366                  priorities=NULL;
4367                  free(branchDirection);
4368                  branchDirection=NULL;
4369                  free(pseudoDown);
4370                  pseudoDown=NULL;
4371                  free(pseudoUp);
4372                  pseudoUp=NULL;
4373                  free(solutionIn);
4374                  solutionIn=NULL;
4375                  free(prioritiesIn);
4376                  prioritiesIn=NULL;
4377                  free(sosStart);
4378                  sosStart=NULL;
4379                  free(sosIndices);
4380                  sosIndices=NULL;
4381                  free(sosType);
4382                  sosType=NULL;
4383                  free(sosReference);
4384                  sosReference=NULL;
4385                  free(cut);
4386                  cut=NULL;
4387                  free(sosPriority);
4388                  sosPriority=NULL;
4389#ifdef COIN_HAS_ASL
4390                }
4391#endif               
4392                if (nodeStrategy) {
4393                  // change default
4394                  if (nodeStrategy>2) {
4395                    // up or down
4396                    int way = (((nodeStrategy-1)%1)==1) ? -1 : +1;
4397                    babModel->setPreferredWay(way);
4398#if 0
4399                    OsiObject ** objects = babModel->objects();
4400                    int numberObjects = babModel->numberObjects();
4401                    for (int iObj = 0;iObj<numberObjects;iObj++) {
4402                      CbcObject * obj =
4403                        dynamic_cast <CbcObject *>(objects[iObj]) ;
4404                      assert (obj);
4405                      obj->setPreferredWay(way);
4406                    }
4407#endif
4408                  }
4409                  if (nodeStrategy==2||nodeStrategy>4) {
4410                    // depth
4411                    CbcCompareDefault compare;
4412                    compare.setWeight(-3.0);
4413                    babModel->setNodeComparison(compare);
4414                  } else if (nodeStrategy==0) {
4415                    // hybrid was default i.e. mixture of low depth and infeasibility
4416                  } else if (nodeStrategy==1) {
4417                    // real fewest
4418                    CbcCompareDefault compare;
4419                    compare.setWeight(-2.0);
4420                    babModel->setNodeComparison(compare);
4421                  }
4422                }
4423                if (cppValue>=0) {
4424                  int prepro = useStrategy ? -1 : preProcess;
4425                  // generate code
4426                  FILE * fp = fopen("user_driver.cpp","w");
4427                  if (fp) {
4428                    // generate enough to do BAB
4429                    babModel->generateCpp(fp,1);
4430                    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4431                    // Make general so do factorization
4432                    int factor = osiclp->getModelPtr()->factorizationFrequency();
4433                    osiclp->getModelPtr()->setFactorizationFrequency(200);
4434                    osiclp->generateCpp(fp);
4435                    osiclp->getModelPtr()->setFactorizationFrequency(factor);
4436                    //solveOptions.generateCpp(fp);
4437                    fclose(fp);
4438                    // now call generate code
4439                    generateCode(babModel,"user_driver.cpp",cppValue,prepro);
4440                  } else {
4441                    std::cout<<"Unable to open file user_driver.cpp"<<std::endl;
4442                  }
4443                }
4444                if (!babModel->numberStrong())
4445                  babModel->setNumberBeforeTrust(0);
4446                if (useStrategy) {
4447                  CbcStrategyDefault strategy(true,babModel->numberStrong(),babModel->numberBeforeTrust());
4448                  strategy.setupPreProcessing(1);
4449                  babModel->setStrategy(strategy);
4450                }
4451                if (testOsiOptions>=0) {
4452                  sprintf(generalPrint,"Testing OsiObject options %d\n",testOsiOptions);
4453                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
4454                    << generalPrint
4455                    <<CoinMessageEol;
4456                  if (!numberSOS) {
4457                    babModel->solver()->findIntegersAndSOS(false);
4458#ifdef COIN_HAS_LINK
4459                    // If linked then pass in model
4460                    OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
4461                    if (solver3) {
4462                      CbcHeuristicDynamic3 serendipity(*babModel);
4463                      serendipity.setHeuristicName("linked");
4464                      babModel->addHeuristic(&serendipity);
4465                      int options = parameters[whichParam(MIPOPTIONS,numberParameters,parameters)].intValue()/10000;
4466                      CglStored stored;
4467                      if (options) {
4468                        printf("nlp options %d\n",options);
4469                        /*
4470                          1 - force mini branch and bound
4471                          2 - set priorities high on continuous
4472                          4 - try adding OA cuts
4473                          8 - try doing quadratic linearization
4474                          16 - try expanding knapsacks
4475                          32 - OA cuts strictly concave
4476                          64 - no branching at all on bilinear x-x!
4477                        */
4478                        if ((options&2)) {
4479                          solver3->setBiLinearPriorities(10,tightenFactor > 0.0 ? tightenFactor : 1.0);
4480                        } else if (tightenFactor>0.0) {
4481                          // set grid size for all continuous bi-linear
4482                          solver3->setMeshSizes(tightenFactor);
4483                        }
4484                        if ((options&4)) {
4485                          solver3->setSpecialOptions2(solver3->specialOptions2()|(8+4));
4486                          // say convex
4487                          solver3->sayConvex((options&32)==0);
4488                        }
4489                        int extra1 = parameters[whichParam(EXTRA1,numberParameters,parameters)].intValue();
4490                        if ((options&1)!=0&&extra1>0)
4491                          solver3->setFixedPriority(extra1);
4492                        double cutoff=COIN_DBL_MAX;
4493                        if ((options&8))
4494                          cutoff=solver3->linearizedBAB(&stored);
4495                        if (cutoff<babModel->getCutoff()) {
4496                          babModel->setCutoff(cutoff);
4497                          // and solution
4498                          //babModel->setBestObjectiveValue(solver3->bestObjectiveValue());
4499                          babModel->setBestSolution(solver3->bestSolution(),solver3->getNumCols(),
4500                                                    solver3->bestObjectiveValue());
4501                        }
4502                        if ((options&64))
4503                          solver3->setBranchingStrategyOnVariables(16,-1,4);
4504                      }
4505                      solver3->setCbcModel(babModel);
4506                      if (stored.sizeRowCuts()) 
4507                        babModel->addCutGenerator(&stored,1,"Stored");
4508                      CglTemporary temp;
4509                      babModel->addCutGenerator(&temp,1,"OnceOnly");
4510                      //choose.setNumberBeforeTrusted(2000);
4511                      //choose.setNumberStrong(20);
4512                    }
4513                    // For temporary testing of heuristics
4514                    //int testOsiOptions = parameters[whichParam(TESTOSI,numberParameters,parameters)].intValue();
4515                    if (testOsiOptions>=10) {
4516                      if (testOsiOptions>=20)
4517                        testOsiOptions -= 10;
4518                      printf("*** Temp heuristic with mode %d\n",testOsiOptions-10);
4519                      OsiSolverLink * solver3 = dynamic_cast<OsiSolverLink *> (babModel->solver());
4520                      assert (solver3) ;
4521                      int extra1 = parameters[whichParam(EXTRA1,numberParameters,parameters)].intValue();
4522                      solver3->setBiLinearPriority(extra1);
4523                      printf("bilinear priority now %d\n",extra1);
4524                      int extra2 = parameters[whichParam(EXTRA2,numberParameters,parameters)].intValue();
4525                      double saveDefault = solver3->defaultBound();
4526                      solver3->setDefaultBound((double) extra2);
4527                      double * solution = solver3->heuristicSolution(slpValue>0 ? slpValue : 40 ,1.0e-5,testOsiOptions-10);
4528                      solver3->setDefaultBound(saveDefault);
4529                      if (!solution)
4530                        printf("Heuristic failed\n");
4531                    }
4532#endif
4533                  } else {
4534                    // move across
4535                    babModel->deleteObjects(false);
4536                    //babModel->addObjects(babModel->solver()->numberObjects(),babModel->solver()->objects());
4537                  }
4538                  CbcBranchDefaultDecision decision;
4539                  if (babModel->numberStrong()) {
4540                    OsiChooseStrong choose(babModel->solver());
4541                    choose.setNumberBeforeTrusted(babModel->numberBeforeTrust());
4542                    choose.setNumberStrong(babModel->numberStrong());
4543                    choose.setShadowPriceMode(testOsiOptions);
4544                    decision.setChooseMethod(choose);
4545                  } else {
4546                    OsiChooseVariable choose(babModel->solver());
4547                    decision.setChooseMethod(choose);
4548                  }
4549                  babModel->setBranchingMethod(decision);
4550                  if (useCosts&&testOsiOptions>=0) {
4551                    int numberColumns = babModel->getNumCols();
4552                    int * sort = new int[numberColumns];
4553                    double * dsort = new double[numberColumns];
4554                    int * priority = new int [numberColumns];
4555                    const double * objective = babModel->getObjCoefficients();
4556                    const double * lower = babModel->getColLower() ;
4557                    const double * upper = babModel->getColUpper() ;
4558                    const CoinPackedMatrix * matrix = babModel->solver()->getMatrixByCol();
4559                    const int * columnLength = matrix->getVectorLengths();
4560                    int iColumn;
4561                    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4562                      sort[iColumn]=iColumn;
4563                      if (useCosts==1)
4564                        dsort[iColumn]=-fabs(objective[iColumn]);
4565                      else if (useCosts==2)
4566                        dsort[iColumn]=iColumn;
4567                      else if (useCosts==3)
4568                        dsort[iColumn]=upper[iColumn]-lower[iColumn];
4569                      else if (useCosts==4)
4570                        dsort[iColumn]=-(upper[iColumn]-lower[iColumn]);
4571                      else if (useCosts==5)
4572                        dsort[iColumn]=-columnLength[iColumn];
4573                    }
4574                    CoinSort_2(dsort,dsort+numberColumns,sort);
4575                    int level=0;
4576                    double last = -1.0e100;
4577                    for (int i=0;i<numberColumns;i++) {
4578                      int iPut=sort[i];
4579                      if (dsort[i]!=last) {
4580                        level++;
4581                        last=dsort[i];
4582                      }
4583                      priority[iPut]=level;
4584                    }
4585                    OsiObject ** objects = babModel->objects();
4586                    int numberObjects = babModel->numberObjects();
4587                    for (int iObj = 0;iObj<numberObjects;iObj++) {
4588                      OsiObject * obj = objects[iObj] ;
4589                      int iColumn = obj->columnNumber();
4590                      if (iColumn>=0)
4591                        obj->setPriority(priority[iColumn]);
4592                    }
4593                    delete [] priority;
4594                    delete [] sort;
4595                    delete [] dsort;
4596                  }
4597                }
4598                checkSOS(babModel, babModel->solver());
4599                if (doSprint>0) {
4600                  // Sprint for primal solves
4601                  ClpSolve::SolveType method = ClpSolve::usePrimalorSprint;
4602                  ClpSolve::PresolveType presolveType = ClpSolve::presolveOff;
4603                  int numberPasses = 5;
4604                  int options[] = {0,3,0,0,0,0};
4605                  int extraInfo[] = {-1,20,-1,-1,-1,-1};
4606                  extraInfo[1]=doSprint;
4607                  int independentOptions[] = {0,0,3};
4608                  ClpSolve clpSolve(method,presolveType,numberPasses,
4609                                    options,extraInfo,independentOptions);
4610                  // say use in OsiClp
4611                  clpSolve.setSpecialOption(6,1);
4612                  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4613                  osiclp->setSolveOptions(clpSolve);
4614                  osiclp->setHintParam(OsiDoDualInResolve,false);
4615                  // switch off row copy
4616                  osiclp->getModelPtr()->setSpecialOptions(osiclp->getModelPtr()->specialOptions()|256);
4617                  osiclp->getModelPtr()->setInfeasibilityCost(1.0e11);
4618                }
4619#ifdef COIN_HAS_LINK
4620                if (storedAmpl.sizeRowCuts()) {
4621                  if (preProcess) {
4622                    const int * originalColumns = process.originalColumns();
4623                    int numberColumns = babModel->getNumCols();
4624                    int * newColumn = new int[numberOriginalColumns];
4625                    int i;
4626                    for (i=0;i<numberOriginalColumns;i++) 
4627                      newColumn[i]=-1;
4628                    for (i=0;i<numberColumns;i++) {
4629                      int iColumn = originalColumns[i];
4630                      newColumn[iColumn]=i;
4631                    }
4632                    int * buildColumn = new int[numberColumns];
4633                    // Build up valid cuts
4634                    int nBad=0;
4635                    int nCuts = storedAmpl.sizeRowCuts();
4636                    CglStored newCuts;
4637                    for (i=0;i<nCuts;i++) {
4638                      const OsiRowCut * cut = storedAmpl.rowCutPointer(i);
4639                      double lb = cut->lb();
4640                      double ub = cut->ub();
4641                      int n=cut->row().getNumElements();
4642                      const int * column = cut->row().getIndices();
4643                      const double * element = cut->row().getElements();
4644                      bool bad=false;
4645                      for (int i=0;i<n;i++) {
4646                        int iColumn = column[i];
4647                        iColumn = newColumn[iColumn];
4648                        if (iColumn>=0) {
4649                          buildColumn[i]=iColumn;
4650                        } else {
4651                          bad=true;
4652                          break;
4653                        }
4654                      }
4655                      if (!bad) {
4656                        newCuts.addCut(lb,ub,n,buildColumn,element);
4657                      } else {
4658                        nBad++;
4659                      }
4660                    }
4661                    storedAmpl=newCuts;
4662                    if (nBad)
4663                      printf("%d cuts dropped\n",nBad);
4664                    delete [] newColumn;
4665                    delete [] buildColumn;
4666                  }
4667                }
4668#endif
4669#ifdef CLP_MALLOC_STATISTICS
4670                malloc_stats();
4671                malloc_stats2();
4672#endif
4673                if (outputFormat==5) {
4674                  osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4675                  lpSolver = osiclp->getModelPtr();
4676                  lpSolver->setPersistenceFlag(1);
4677                }
4678#ifdef COIN_HAS_ASL
4679                // add in lotsizing
4680                if (usingAmpl&&info.special) {
4681                  int numberColumns = babModel->getNumCols();
4682                  int i;
4683                  int n=0;
4684                  if (preProcess) {
4685                    const int * originalColumns = process.originalColumns();
4686                    for (i=0;i<numberColumns;i++) {
4687                      int iColumn = originalColumns[i];
4688                      assert (iColumn>=i);
4689                      int iType = info.special[iColumn];
4690                      if (iType) {
4691                        assert (iType==1);
4692                        n++;
4693                      }
4694                      info.special[i]=iType;
4695                    }
4696                  }
4697                  if (n) {
4698                    int numberIntegers=0;
4699                    int numberOldObjects=0;
4700                    OsiObject ** oldObjects=NULL;
4701                    const double * lower = babModel->solver()->getColLower();
4702                    const double * upper = babModel->solver()->getColUpper();
4703                    if (testOsiOptions<0) {
4704                      // *************************************************************
4705                      // CbcObjects
4706                      numberIntegers = babModel->numberIntegers();
4707                      /* model may not have created objects
4708                         If none then create
4709                      */
4710                      if (!numberIntegers||!babModel->numberObjects()) {
4711                        int type = (pseudoUp) ? 1 : 0;
4712                        babModel->findIntegers(true,type);
4713                        numberIntegers = babModel->numberIntegers();
4714                      }
4715                      oldObjects = babModel->objects();
4716                      numberOldObjects = babModel->numberObjects();
4717                    } else {
4718                      numberIntegers = testOsiSolver->getNumIntegers();
4719                      if (!numberIntegers||!testOsiSolver->numberObjects()) {
4720                        /* model may not have created objects
4721                           If none then create
4722                        */
4723                        testOsiSolver->findIntegers(false);
4724                        numberIntegers = testOsiSolver->getNumIntegers();
4725                      }
4726                      oldObjects = testOsiSolver->objects();
4727                      numberOldObjects = testOsiSolver->numberObjects();
4728                    }
4729                    OsiObject ** objects = new OsiObject * [n];
4730                    n=0;
4731                    // set new objects to have one lower priority
4732                    double ranges[] = {-COIN_DBL_MAX,-1.0,1.0,COIN_DBL_MAX};
4733                    for (int iObj = 0;iObj<numberOldObjects;iObj++) {
4734                      int iColumn = oldObjects[iObj]->columnNumber();
4735                      if (iColumn>=0&&info.special[iColumn]) {
4736                        if (lower[iColumn]<=-1.0&&upper[iColumn]>=0.0) {
4737                          ranges[0]=lower[iColumn];
4738                          ranges[3]=upper[iColumn];
4739                          int priority = oldObjects[iObj]->priority();
4740                          if (testOsiOptions<0) {
4741                            objects[n] = new CbcLotsize(babModel,iColumn,2,ranges,true);
4742                          } else {
4743                            objects[n] = new OsiLotsize(testOsiSolver,iColumn,2,ranges,true);
4744                          }
4745                          objects[n++]->setPriority (priority-1);
4746                        }
4747                      }
4748                    }
4749                    if (testOsiOptions<0) {
4750                      babModel->addObjects(n,objects);
4751                    } else {
4752                      testOsiSolver->addObjects(n,objects);
4753                    }
4754                    for (i=0;i<n;i++)
4755                      delete objects[i];
4756                    delete [] objects;
4757                  }
4758                }
4759#endif
4760                if (storedAmpl.sizeRowCuts()) {
4761                  //babModel->addCutGenerator(&storedAmpl,1,"AmplStored");
4762                  int numberRowCuts = storedAmpl.sizeRowCuts();
4763                  for (int i=0;i<numberRowCuts;i++) {
4764                    const OsiRowCut * rowCutPointer = storedAmpl.rowCutPointer(i);
4765                    babModel->makeGlobalCut(rowCutPointer);
4766                  }
4767                }
4768                // If defaults then increase trust for small models
4769                if (!strongChanged) {
4770                  int numberColumns = babModel->getNumCols();
4771                  if (numberColumns<=50)
4772                    babModel->setNumberBeforeTrust(1000);
4773                  else if (numberColumns<=100)
4774                    babModel->setNumberBeforeTrust(100);
4775                  else if (numberColumns<=300)
4776                    babModel->setNumberBeforeTrust(50);
4777                }
4778#ifdef CBC_THREAD
4779                int numberThreads =parameters[whichParam(THREADS,numberParameters,parameters)].intValue();
4780                babModel->setNumberThreads(numberThreads%100);
4781                babModel->setThreadMode(numberThreads/100);
4782#endif
4783                int returnCode=callBack(babModel,3);
4784                if (returnCode) {
4785                  // exit if user wants
4786                  delete babModel;
4787                  return returnCode;
4788                }
4789                babModel->branchAndBound(statistics);
4790                returnCode=callBack(babModel,4);
4791                if (returnCode) {
4792                  // exit if user wants
4793                  model.moveInfo(*babModel);
4794                  delete babModel;
4795                  return returnCode;
4796                }
4797#ifdef CLP_MALLOC_STATISTICS
4798                malloc_stats();
4799                malloc_stats2();
4800#endif
4801                checkSOS(babModel, babModel->solver());
4802              } else if (type==MIPLIB) {
4803                CbcStrategyDefault strategy(true,babModel->numberStrong(),babModel->numberBeforeTrust());
4804                // Set up pre-processing
4805                int translate2[]={9999,1,1,3,2,4,5};
4806                if (preProcess)
4807                  strategy.setupPreProcessing(translate2[preProcess]);
4808                babModel->setStrategy(strategy);
4809#ifdef CBC_THREAD
4810                int numberThreads =parameters[whichParam(THREADS,numberParameters,parameters)].intValue();
4811                babModel->setNumberThreads(numberThreads%100);
4812                babModel->setThreadMode(numberThreads/100);
4813#endif
4814                if (outputFormat==5) {
4815                  osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4816                  lpSolver = osiclp->getModelPtr();
4817                  lpSolver->setPersistenceFlag(1);
4818                }
4819                if (testOsiOptions>=0) {
4820                  printf("Testing OsiObject options %d\n",testOsiOptions);
4821                  CbcBranchDefaultDecision decision;
4822                  OsiChooseStrong choose(babModel->solver());
4823                  choose.setNumberBeforeTrusted(babModel->numberBeforeTrust());
4824                  choose.setNumberStrong(babModel->numberStrong());
4825                  choose.setShadowPriceMode(testOsiOptions);
4826                  //babModel->deleteObjects(false);
4827                  decision.setChooseMethod(choose);
4828                  babModel->setBranchingMethod(decision);
4829                }
4830                model = *babModel;
4831                return 777;
4832              } else {
4833                strengthenedModel = babModel->strengthenedModel();
4834              }
4835              currentBranchModel = NULL;
4836              osiclp = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4837              if (debugFile=="createAfterPre"&&babModel->bestSolution()) {
4838                lpSolver = osiclp->getModelPtr();
4839                //move best solution (should be there -- but ..)
4840                int n = lpSolver->getNumCols();
4841                memcpy(lpSolver->primalColumnSolution(),babModel->bestSolution(),n*sizeof(double));
4842                saveSolution(osiclp->getModelPtr(),"debug.file");
4843              }
4844              if (!noPrinting) {
4845                // Print more statistics
4846                sprintf(generalPrint,"Cuts at root node changed objective from %g to %g",
4847                        babModel->getContinuousObjective(),babModel->rootObjectiveAfterCuts());
4848                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4849                  << generalPrint
4850                  <<CoinMessageEol;
4851               
4852                numberGenerators = babModel->numberCutGenerators();
4853                char timing[30];
4854                for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
4855                  CbcCutGenerator * generator = babModel->cutGenerator(iGenerator);
4856                  sprintf(generalPrint,"%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts",
4857                          generator->cutGeneratorName(),
4858                          generator->numberTimesEntered(),
4859                          generator->numberCutsInTotal(),
4860                          generator->numberCutsActive());
4861                  if (generator->timing()) {
4862                    sprintf(timing," (%.3f seconds)",generator->timeInCutGenerator());
4863                    strcat(generalPrint,timing);
4864                  }
4865                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
4866                    << generalPrint
4867                    <<CoinMessageEol;
4868                }
4869              }
4870              // adjust time to allow for children on some systems
4871              time2 = CoinCpuTime() + CoinCpuTimeJustChildren();
4872              totalTime += time2-time1;
4873              // For best solution
4874              double * bestSolution = NULL;
4875              if (babModel->getMinimizationObjValue()<1.0e50&&type==BAB) {
4876                // post process
4877                int n;
4878                if (preProcess) {
4879                  n = saveSolver->getNumCols();
4880                  bestSolution = new double [n];
4881                  OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4882                  ClpSimplex * lpSolver = clpSolver->getModelPtr();
4883                  lpSolver->setSpecialOptions(lpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
4884                  process.postProcess(*babModel->solver());
4885                  // Solution now back in saveSolver
4886                  babModel->assignSolver(saveSolver);
4887                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
4888                } else {
4889                  n = babModel->solver()->getNumCols();
4890                  bestSolution = new double [n];
4891                  memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
4892                }
4893                model.setBestSolution(bestSolution,n,babModel->getObjValue());
4894                // and put back in very original solver
4895                {
4896                  ClpSimplex * original = originalSolver->getModelPtr();
4897                  double * lower = original->columnLower();
4898                  double * upper = original->columnUpper();
4899                  double * solution = original->primalColumnSolution();
4900                  int n = original->numberColumns();
4901                  //assert (!n||n==babModel->solver()->getNumCols());
4902                  for (int i=0;i<n;i++) {
4903                    solution[i]=bestSolution[i];
4904                    if (originalSolver->isInteger(i)) {
4905                      lower[i]=solution[i];
4906                      upper[i]=solution[i];
4907                    }
4908                  }
4909                }
4910                checkSOS(babModel, babModel->solver());
4911              }
4912              if (type==STRENGTHEN&&strengthenedModel)
4913                clpSolver = dynamic_cast< OsiClpSolverInterface*> (strengthenedModel);
4914#ifdef COIN_HAS_ASL
4915              else if (usingAmpl) 
4916                clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4917#endif
4918              lpSolver = clpSolver->getModelPtr();
4919              if (numberChanged) {
4920                for (int i=0;i<numberChanged;i++) {
4921                  int iColumn=changed[i];
4922                  clpSolver->setContinuous(iColumn);
4923                }
4924                delete [] changed;
4925              }
4926              if (type==BAB) {
4927                //move best solution (should be there -- but ..)
4928                int n = lpSolver->getNumCols();
4929                if (bestSolution) {
4930                  memcpy(lpSolver->primalColumnSolution(),bestSolution,n*sizeof(double));
4931                  // now see what that does to row solution
4932                  int numberRows=lpSolver->numberRows();
4933                  double * rowSolution = lpSolver->primalRowSolution();
4934                  memset (rowSolution,0,numberRows*sizeof(double));
4935                  lpSolver->clpMatrix()->times(1.0,bestSolution,rowSolution);
4936                  lpSolver->setObjectiveValue(babModel->getObjValue());
4937                }
4938                if (debugFile=="create"&&bestSolution) {
4939                  saveSolution(lpSolver,"debug.file");
4940                }
4941                delete [] bestSolution;
4942                std::string statusName[]={"Finished","Stopped on ","Difficulties",
4943                                          "","","User ctrl-c"};
4944                std::string minor[]={"","","gap","nodes","time","","solutions","user ctrl-c"};
4945                int iStat = babModel->status();
4946                int iStat2 = babModel->secondaryStatus();
4947                if (!noPrinting) {
4948                  sprintf(generalPrint,"Result - %s%s objective %.16g after %d nodes and %d iterations - took %.2f seconds",
4949                          statusName[iStat].c_str(),minor[iStat2].c_str(),
4950                          babModel->getObjValue(),babModel->getNodeCount(),
4951                          babModel->getIterationCount(),time2-time1);
4952                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
4953                    << generalPrint
4954                    <<CoinMessageEol;
4955                }
4956                int returnCode=callBack(babModel,5);
4957                if (returnCode) {
4958                  // exit if user wants
4959                  model.moveInfo(*babModel);
4960                  delete babModel;
4961                  return returnCode;
4962                }
4963#ifdef COIN_HAS_ASL
4964                if (usingAmpl) {
4965                  clpSolver = dynamic_cast< OsiClpSolverInterface*> (babModel->solver());
4966                  lpSolver = clpSolver->getModelPtr();
4967                  double value = babModel->getObjValue()*lpSolver->getObjSense();
4968                  char buf[300];
4969                  int pos=0;
4970                  if (iStat==0) {
4971                    if (babModel->getObjValue()<1.0e40) {
4972                      pos += sprintf(buf+pos,"optimal," );
4973                    } else {
4974                      // infeasible
4975                      iStat=1;
4976                      pos += sprintf(buf+pos,"infeasible,");
4977                    }
4978                  } else if (iStat==1) {
4979                    if (iStat2!=6)
4980                      iStat=3;
4981                    else
4982                      iStat=4;
4983                    pos += sprintf(buf+pos,"stopped on %s,",minor[iStat2].c_str());
4984                  } else if (iStat==2) {
4985                    iStat = 7;
4986                    pos += sprintf(buf+pos,"stopped on difficulties,");
4987                  } else if (iStat==5) {
4988                    iStat = 3;
4989                    pos += sprintf(buf+pos,"stopped on ctrl-c,");
4990                  } else {
4991                    pos += sprintf(buf+pos,"status unknown,");
4992                    iStat=6;
4993                  }
4994                  info.problemStatus=iStat;
4995                  info.objValue = value;
4996                  if (babModel->getObjValue()<1.0e40) {
4997                    int precision = ampl_obj_prec();
4998                    if (precision>0)
4999                      pos += sprintf(buf+pos," objective %.*g",precision,
5000                                     value);
5001                    else
5002                      pos += sprintf(buf+pos," objective %g",value);
5003                  }
5004                  sprintf(buf+pos,"\n%d nodes, %d iterations, %g seconds",
5005                          babModel->getNodeCount(),
5006                          babModel->getIterationCount(),
5007                          totalTime);
5008                  if (bestSolution) {
5009                    free(info.primalSolution);
5010                    if (!numberKnapsack) {
5011                      info.primalSolution = (double *) malloc(n*sizeof(double));
5012                      CoinCopyN(lpSolver->primalColumnSolution(),n,info.primalSolution);
5013                      int numberRows = lpSolver->numberRows();
5014                      free(info.dualSolution);
5015                      info.dualSolution = (double *) malloc(numberRows*sizeof(double));
5016                      CoinCopyN(lpSolver->dualRowSolution(),numberRows,info.dualSolution);
5017                    } else {
5018                      // expanded knapsack
5019                      info.dualSolution=NULL;
5020                      int numberColumns = saveCoinModel.numberColumns();
5021                      info.primalSolution = (double *) malloc(numberColumns*sizeof(double));
5022                      // Fills in original solution (coinModel length)
5023                      afterKnapsack(saveTightenedModel,  whichColumn,  knapsackStart, 
5024                                    knapsackRow,  numberKnapsack,
5025                                    lpSolver->primalColumnSolution(), info.primalSolution,1);
5026                    }
5027                  } else {
5028                    info.primalSolution=NULL;
5029                    info.dualSolution=NULL;
5030                  }
5031                  // put buffer into info
5032                  strcpy(info.buffer,buf);
5033                }
5034#endif
5035              } else {
5036                std::cout<<"Model strengthened - now has "<<clpSolver->getNumRows()
5037                         <<" rows"<<std::endl;
5038              }
5039              time1 = time2;