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

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

fpump - make dj fixing optional

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