source: branches/BSP/trunk/Cbc/src/CbcSolver.cpp @ 762

Last change on this file since 762 was 762, checked in by tkr, 14 years ago

Fixing dependence of libCbc on libCbcSolver

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