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

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

comment fake stuff

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