source: branches/devel/Cbc/src/CoinSolve.cpp @ 578

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

for Jon

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