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

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

comment out debug malloc!

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