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

Last change on this file since 759 was 759, checked in by forrest, 14 years ago

start of work on new vub heuristic

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