source: trunk/Cbc/src/CbcSolver.cpp @ 857

Last change on this file since 857 was 857, checked in by forrest, 11 years ago

allow dense solver

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