source: stable/2.0/Cbc/src/CbcSolver.cpp @ 844

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

make useSolution work better

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