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

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

keep up to date

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