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

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

try and speed up feasibility pump

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