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

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

save basis in solver correctly

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