source: branches/devel/Cbc/src/CbcSolver.cpp @ 612

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

for nonlinear

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