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

Last change on this file since 763 was 763, checked in by andreasw, 14 years ago

unded [762] and moved CbcCbcParam?.cpp from libCbc to libCbcSolver

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