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

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

messages for Stefan

File size: 326.2 KB
Line 
1// Copyright (C) 2007, International Business Machines
2// Corporation and others.  All Rights Reserved.
3   
4#include "CbcConfig.h"
5#include "CoinPragma.hpp"
6
7#include <cassert>
8#include <cstdio>
9#include <cmath>
10#include <cfloat>
11#include <cstring>
12#include <iostream>
13
14
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              int returnCode;
4026              const char * message = 
4027                parameters_[iParam].setDoubleParameterWithMessage(lpSolver,value,returnCode);
4028              if (!noPrinting_&&strlen(message)) {
4029                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4030                  << message
4031                  <<CoinMessageEol;
4032              }
4033            } else if (type<81) {
4034              int returnCode;
4035              const char * message = 
4036                parameters_[iParam].setDoubleParameterWithMessage(model_,value,returnCode);
4037              if (!noPrinting_&&strlen(message)) {
4038                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4039                  << message
4040                  <<CoinMessageEol;
4041              }
4042            } else {
4043              int returnCode;
4044              const char * message = 
4045                parameters_[iParam].setDoubleParameterWithMessage(lpSolver,value,returnCode);
4046              if (!noPrinting_&&strlen(message)) {
4047                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4048                  << message
4049                  <<CoinMessageEol;
4050              }
4051              switch(type) {
4052              case DJFIX:
4053                djFix=value;
4054                if (goodModel&&djFix<1.0e20) {
4055                  // do some fixing
4056                  clpSolver = dynamic_cast< OsiClpSolverInterface*> (model_.solver());
4057                  clpSolver->initialSolve();
4058                  lpSolver = clpSolver->getModelPtr();
4059                  int numberColumns = lpSolver->numberColumns();
4060                  int i;
4061                  const char * type = lpSolver->integerInformation();
4062                  double * lower = lpSolver->columnLower();
4063                  double * upper = lpSolver->columnUpper();
4064                  double * solution = lpSolver->primalColumnSolution();
4065                  double * dj = lpSolver->dualColumnSolution();
4066                  int numberFixed=0;
4067                  double dextra4 = parameters_[whichParam(DEXTRA4,numberParameters_,parameters_)].doubleValue();
4068                  if (dextra4)
4069                    printf("Multiple for continuous dj fixing is %g\n",dextra4);
4070                  for (i=0;i<numberColumns;i++) {
4071                    double djValue = dj[i];
4072                    if (!type[i])
4073                      djValue *= dextra4;
4074                    if (type[i]||dextra4) {
4075                      double value = solution[i];
4076                      if (value<lower[i]+1.0e-5&&djValue>djFix) {
4077                        solution[i]=lower[i];
4078                        upper[i]=lower[i];
4079                        numberFixed++;
4080                      } else if (value>upper[i]-1.0e-5&&djValue<-djFix) {
4081                        solution[i]=upper[i];
4082                        lower[i]=upper[i];
4083                        numberFixed++;
4084                      }
4085                    }
4086                  }
4087                  sprintf(generalPrint,"%d columns fixed\n",numberFixed);
4088                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
4089                    << generalPrint
4090                    <<CoinMessageEol;
4091                }
4092                break;
4093              case GAPRATIO:
4094                gapRatio=value;
4095                break;
4096              case TIGHTENFACTOR:
4097                tightenFactor=value;
4098                if(!complicatedInteger)
4099                  defaultSettings=false; // user knows what she is doing
4100                break;
4101              default:
4102                break;
4103              }
4104            }
4105          } else if (valid==1) {
4106            std::cout<<" is illegal for double parameter "<<parameters_[iParam].name()<<" value remains "<<
4107              parameters_[iParam].doubleValue()<<std::endl;
4108          } else {
4109            std::cout<<parameters_[iParam].name()<<" has value "<<
4110              parameters_[iParam].doubleValue()<<std::endl;
4111          }
4112        } else if (type<201) {
4113          // get next field as int
4114          int value = CoinReadGetIntField(argc,argv,&valid);
4115          if (!valid) {
4116            if (type<151) {
4117              if (parameters_[iParam].type()==PRESOLVEPASS)
4118                preSolve = value;
4119              else if (parameters_[iParam].type()==IDIOT)
4120                doIdiot = value;
4121              else if (parameters_[iParam].type()==SPRINT)
4122                doSprint = value;
4123              else if (parameters_[iParam].type()==OUTPUTFORMAT)
4124                outputFormat = value;
4125              else if (parameters_[iParam].type()==SLPVALUE)
4126                slpValue = value;
4127              else if (parameters_[iParam].type()==CPP)
4128                cppValue = value;
4129              else if (parameters_[iParam].type()==PRESOLVEOPTIONS)
4130                presolveOptions = value;
4131              else if (parameters_[iParam].type()==PRINTOPTIONS)
4132                printOptions = value;
4133              else if (parameters_[iParam].type()==SUBSTITUTION)
4134                substitution = value;
4135              else if (parameters_[iParam].type()==DUALIZE)
4136                dualize = value;
4137              else if (parameters_[iParam].type()==PROCESSTUNE)
4138                tunePreProcess = value;
4139              else if (parameters_[iParam].type()==USESOLUTION)
4140                useSolution = value;
4141              else if (parameters_[iParam].type()==VERBOSE)
4142                verbose = value;
4143              int returnCode;
4144              const char * message = 
4145                parameters_[iParam].setIntParameterWithMessage(lpSolver,value,returnCode);
4146              if (!noPrinting_&&strlen(message)) {
4147                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4148                  << message
4149                  <<CoinMessageEol;
4150              }
4151            } else {
4152              if (parameters_[iParam].type()==CUTPASS)
4153                cutPass = value;
4154              else if (parameters_[iParam].type()==CUTPASSINTREE)
4155                cutPassInTree = value;
4156              else if (parameters_[iParam].type()==STRONGBRANCHING||
4157                       parameters_[iParam].type()==NUMBERBEFORE)
4158                strongChanged=true;
4159              int returnCode;
4160              const char * message = 
4161                parameters_[iParam].setIntParameterWithMessage(model_,value,returnCode);
4162              if (!noPrinting_&&strlen(message)) {
4163                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4164                  << message
4165                  <<CoinMessageEol;
4166              }
4167            }
4168          } else if (valid==1) {
4169            std::cout<<" is illegal for integer parameter "<<parameters_[iParam].name()<<" value remains "<<
4170              parameters_[iParam].intValue()<<std::endl;
4171          } else {
4172            std::cout<<parameters_[iParam].name()<<" has value "<<
4173              parameters_[iParam].intValue()<<std::endl;
4174          }
4175        } else if (type<301) {
4176          // one of several strings
4177          std::string value = CoinReadGetString(argc,argv);
4178          int action = parameters_[iParam].parameterOption(value);
4179          if (action<0) {
4180            if (value!="EOL") {
4181              // no match
4182              parameters_[iParam].printOptions();
4183            } else {
4184              // print current value
4185              std::cout<<parameters_[iParam].name()<<" has value "<<
4186                parameters_[iParam].currentOption()<<std::endl;
4187            }
4188          } else {
4189            const char * message = 
4190              parameters_[iParam].setCurrentOptionWithMessage(action);
4191            if (!noPrinting_) {
4192              generalMessageHandler->message(CLP_GENERAL,generalMessages)
4193                << message
4194                <<CoinMessageEol;
4195            }
4196            // for now hard wired
4197            switch (type) {
4198            case DIRECTION:
4199              if (action==0)
4200                lpSolver->setOptimizationDirection(1);
4201              else if (action==1)
4202                lpSolver->setOptimizationDirection(-1);
4203              else
4204                lpSolver->setOptimizationDirection(0);
4205              break;
4206            case DUALPIVOT:
4207              if (action==0) {
4208                ClpDualRowSteepest steep(3);
4209                lpSolver->setDualRowPivotAlgorithm(steep);
4210              } else if (action==1) {
4211                //ClpDualRowDantzig dantzig;
4212                ClpDualRowSteepest dantzig(5);
4213                lpSolver->setDualRowPivotAlgorithm(dantzig);
4214              } else if (action==2) {
4215                // partial steep
4216                ClpDualRowSteepest steep(2);
4217                lpSolver->setDualRowPivotAlgorithm(steep);
4218              } else {
4219                ClpDualRowSteepest steep;
4220                lpSolver->setDualRowPivotAlgorithm(steep);
4221              }
4222              break;
4223            case PRIMALPIVOT:
4224              if (action==0) {
4225                ClpPrimalColumnSteepest steep(3);
4226                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4227              } else if (action==1) {
4228                ClpPrimalColumnSteepest steep(0);
4229                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4230              } else if (action==2) {
4231                ClpPrimalColumnDantzig dantzig;
4232                lpSolver->setPrimalColumnPivotAlgorithm(dantzig);
4233              } else if (action==3) {
4234                ClpPrimalColumnSteepest steep(4);
4235                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4236              } else if (action==4) {
4237                ClpPrimalColumnSteepest steep(1);
4238                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4239              } else if (action==5) {
4240                ClpPrimalColumnSteepest steep(2);
4241                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4242              } else if (action==6) {
4243                ClpPrimalColumnSteepest steep(10);
4244                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4245              }
4246              break;
4247            case SCALING:
4248              lpSolver->scaling(action);
4249              solver->setHintParam(OsiDoScale,action!=0,OsiHintTry);
4250              doScaling = 1-action;
4251              break;
4252            case AUTOSCALE:
4253              lpSolver->setAutomaticScaling(action!=0);
4254              break;
4255            case SPARSEFACTOR:
4256              lpSolver->setSparseFactorization((1-action)!=0);
4257              break;
4258            case BIASLU:
4259              lpSolver->factorization()->setBiasLU(action);
4260              break;
4261            case PERTURBATION:
4262              if (action==0)
4263                lpSolver->setPerturbation(50);
4264              else
4265                lpSolver->setPerturbation(100);
4266              break;
4267            case ERRORSALLOWED:
4268              allowImportErrors = action;
4269              break;
4270            case INTPRINT:
4271              printMode=action;
4272              break;
4273              //case ALGORITHM:
4274              //algorithm  = action;
4275              //defaultSettings=false; // user knows what she is doing
4276              //abort();
4277              //break;
4278            case KEEPNAMES:
4279              keepImportNames = 1-action;
4280              break;
4281            case PRESOLVE:
4282              if (action==0)
4283                preSolve = 5;
4284              else if (action==1)
4285                preSolve=0;
4286              else if (action==2)
4287                preSolve=10;
4288              else
4289                preSolveFile=true;
4290              break;
4291            case PFI:
4292              lpSolver->factorization()->setForrestTomlin(action==0);
4293              break;
4294            case CRASH:
4295              doCrash=action;
4296              break;
4297            case VECTOR:
4298              doVector=action;
4299              break;
4300            case MESSAGES:
4301              lpSolver->messageHandler()->setPrefix(action!=0);
4302              break;
4303            case CHOLESKY:
4304              choleskyType = action;
4305              break;
4306            case GAMMA:
4307              gamma=action;
4308              break;
4309            case BARRIERSCALE:
4310              scaleBarrier=action;
4311              break;
4312            case KKT:
4313              doKKT=action;
4314              break;
4315            case CROSSOVER:
4316              crossover=action;
4317              break;
4318            case SOS:
4319              doSOS=action;
4320              break;
4321            case GOMORYCUTS:
4322              defaultSettings=false; // user knows what she is doing
4323              gomoryAction = action;
4324              break;
4325            case PROBINGCUTS:
4326              defaultSettings=false; // user knows what she is doing
4327              probingAction = action;
4328              break;
4329            case KNAPSACKCUTS:
4330              defaultSettings=false; // user knows what she is doing
4331              knapsackAction = action;
4332              break;
4333            case REDSPLITCUTS:
4334              defaultSettings=false; // user knows what she is doing
4335              redsplitAction = action;
4336              break;
4337            case CLIQUECUTS:
4338              defaultSettings=false; // user knows what she is doing
4339              cliqueAction = action;
4340              break;
4341            case FLOWCUTS:
4342              defaultSettings=false; // user knows what she is doing
4343              flowAction = action;
4344              break;
4345            case MIXEDCUTS:
4346              defaultSettings=false; // user knows what she is doing
4347              mixedAction = action;
4348              break;
4349            case TWOMIRCUTS:
4350              defaultSettings=false; // user knows what she is doing
4351              twomirAction = action;
4352              break;
4353            case LANDPCUTS:
4354              defaultSettings=false; // user knows what she is doing
4355              landpAction = action;
4356              break;
4357            case RESIDCUTS:
4358              defaultSettings=false; // user knows what she is doing
4359              residualCapacityAction = action;
4360              break;
4361            case ROUNDING:
4362              defaultSettings=false; // user knows what she is doing
4363              break;
4364            case FPUMP:
4365              defaultSettings=false; // user knows what she is doing
4366              break;
4367            case RINS:
4368              break;
4369            case CUTSSTRATEGY:
4370              gomoryAction = action;
4371              probingAction = action;
4372              knapsackAction = action;
4373              cliqueAction = action;
4374              flowAction = action;
4375              mixedAction = action;
4376              twomirAction = action;
4377              //landpAction = action;
4378              parameters_[whichParam(GOMORYCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4379              parameters_[whichParam(PROBINGCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4380              parameters_[whichParam(KNAPSACKCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4381              parameters_[whichParam(CLIQUECUTS,numberParameters_,parameters_)].setCurrentOption(action);
4382              parameters_[whichParam(FLOWCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4383              parameters_[whichParam(MIXEDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4384              parameters_[whichParam(TWOMIRCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4385              if (!action) {
4386                redsplitAction = action;
4387                parameters_[whichParam(REDSPLITCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4388                landpAction = action;
4389                parameters_[whichParam(LANDPCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4390                residualCapacityAction = action;
4391                parameters_[whichParam(RESIDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4392              }
4393              break;
4394            case HEURISTICSTRATEGY:
4395              parameters_[whichParam(ROUNDING,numberParameters_,parameters_)].setCurrentOption(action);
4396              parameters_[whichParam(GREEDY,numberParameters_,parameters_)].setCurrentOption(action);
4397              parameters_[whichParam(COMBINE,numberParameters_,parameters_)].setCurrentOption(action);
4398              //parameters_[whichParam(LOCALTREE,numberParameters_,parameters_)].setCurrentOption(action);
4399              parameters_[whichParam(FPUMP,numberParameters_,parameters_)].setCurrentOption(action);
4400              break;
4401            case GREEDY:
4402              defaultSettings=false; // user knows what she is doing
4403              break;
4404            case COMBINE:
4405              defaultSettings=false; // user knows what she is doing
4406              break;
4407            case LOCALTREE:
4408              defaultSettings=false; // user knows what she is doing
4409              break;
4410            case COSTSTRATEGY:
4411              useCosts=action;
4412              break;
4413            case NODESTRATEGY:
4414              nodeStrategy=action;
4415              break;
4416            case PREPROCESS:
4417              preProcess = action;
4418              break;
4419            default:
4420              abort();
4421            }
4422          }
4423        } else {
4424          // action
4425          if (type==EXIT) {
4426#ifdef COIN_HAS_ASL
4427            if(statusUserFunction_[0]) {
4428              if (info.numberIntegers||info.numberBinary) {
4429                // integer
4430              } else {
4431                // linear
4432              }
4433              writeAmpl(&info);
4434              freeArrays2(&info);
4435              freeArgs(&info);
4436            }
4437#endif
4438            break; // stop all
4439          }
4440          switch (type) {
4441          case DUALSIMPLEX:
4442          case PRIMALSIMPLEX:
4443          case SOLVECONTINUOUS:
4444          case BARRIER:
4445            if (goodModel) {
4446              double objScale = 
4447                parameters_[whichParam(OBJSCALE2,numberParameters_,parameters_)].doubleValue();
4448              if (objScale!=1.0) {
4449                int iColumn;
4450                int numberColumns=lpSolver->numberColumns();
4451                double * dualColumnSolution = 
4452                  lpSolver->dualColumnSolution();
4453                ClpObjective * obj = lpSolver->objectiveAsObject();
4454                assert(dynamic_cast<ClpLinearObjective *> (obj));
4455                double offset;
4456                double * objective = obj->gradient(NULL,NULL,offset,true);
4457                for (iColumn=0;iColumn<numberColumns;iColumn++) {
4458                  dualColumnSolution[iColumn] *= objScale;
4459                  objective[iColumn] *= objScale;;
4460                }
4461                int iRow;
4462                int numberRows=lpSolver->numberRows();
4463                double * dualRowSolution = 
4464                  lpSolver->dualRowSolution();
4465                for (iRow=0;iRow<numberRows;iRow++) 
4466                  dualRowSolution[iRow] *= objScale;
4467                lpSolver->setObjectiveOffset(objScale*lpSolver->objectiveOffset());
4468              }
4469              ClpSolve::SolveType method;
4470              ClpSolve::PresolveType presolveType;
4471              ClpSimplex * model2 = lpSolver;
4472              if (dualize) {
4473                bool tryIt=true;
4474                double fractionColumn=1.0;
4475                double fractionRow=1.0;
4476                if (dualize==3) {
4477                  dualize=1;
4478                  int numberColumns=lpSolver->numberColumns();
4479                  int numberRows=lpSolver->numberRows();
4480                  if (numberRows<50000||5*numberColumns>numberRows) {
4481                    tryIt=false;
4482                  } else {
4483                    fractionColumn=0.1;
4484                    fractionRow=0.1;
4485                  }
4486                }
4487                if (tryIt) {
4488                  model2 = ((ClpSimplexOther *) model2)->dualOfModel(fractionRow,fractionColumn);
4489                  if (model2) {
4490                    sprintf(generalPrint,"Dual of model has %d rows and %d columns",
4491                            model2->numberRows(),model2->numberColumns());
4492                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4493                      << generalPrint
4494                      <<CoinMessageEol;
4495                    model2->setOptimizationDirection(1.0);
4496                  } else {
4497                    model2 = lpSolver;
4498                    dualize=0;
4499                  }
4500                } else {
4501                  dualize=0;
4502                }
4503              }
4504              if (noPrinting_)
4505                lpSolver->setLogLevel(0);
4506              ClpSolve solveOptions;
4507              solveOptions.setPresolveActions(presolveOptions);
4508              solveOptions.setSubstitution(substitution);
4509              if (preSolve!=5&&preSolve) {
4510                presolveType=ClpSolve::presolveNumber;
4511                if (preSolve<0) {
4512                  preSolve = - preSolve;
4513                  if (preSolve<=100) {
4514                    presolveType=ClpSolve::presolveNumber;
4515                    sprintf(generalPrint,"Doing %d presolve passes - picking up non-costed slacks",
4516                           preSolve);
4517                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4518                      << generalPrint
4519                      <<CoinMessageEol;
4520                    solveOptions.setDoSingletonColumn(true);
4521                  } else {
4522                    preSolve -=100;
4523                    presolveType=ClpSolve::presolveNumberCost;
4524                    sprintf(generalPrint,"Doing %d presolve passes - picking up costed slacks",
4525                           preSolve);
4526                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4527                      << generalPrint
4528                      <<CoinMessageEol;
4529                  }
4530                } 
4531              } else if (preSolve) {
4532                presolveType=ClpSolve::presolveOn;
4533              } else {
4534                presolveType=ClpSolve::presolveOff;
4535              }
4536              solveOptions.setPresolveType(presolveType,preSolve);
4537              if (type==DUALSIMPLEX||type==SOLVECONTINUOUS) {
4538                method=ClpSolve::useDual;
4539              } else if (type==PRIMALSIMPLEX) {
4540                method=ClpSolve::usePrimalorSprint;
4541              } else {
4542                method = ClpSolve::useBarrier;
4543                if (crossover==1) {
4544                  method=ClpSolve::useBarrierNoCross;
4545                } else if (crossover==2) {
4546                  ClpObjective * obj = lpSolver->objectiveAsObject();
4547                  if (obj->type()>1) {
4548                    method=ClpSolve::useBarrierNoCross;
4549                    presolveType=ClpSolve::presolveOff;
4550                    solveOptions.setPresolveType(presolveType,preSolve);
4551                  } 
4552                }
4553              }
4554              solveOptions.setSolveType(method);
4555              if(preSolveFile)
4556                presolveOptions |= 0x40000000;
4557              solveOptions.setSpecialOption(4,presolveOptions);
4558              solveOptions.setSpecialOption(5,printOptions);
4559              if (doVector) {
4560                ClpMatrixBase * matrix = lpSolver->clpMatrix();
4561                if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
4562                  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
4563                  clpMatrix->makeSpecialColumnCopy();
4564                }
4565              }
4566              if (method==ClpSolve::useDual) {
4567                // dual
4568                if (doCrash)
4569                  solveOptions.setSpecialOption(0,1,doCrash); // crash
4570                else if (doIdiot)
4571                  solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
4572              } else if (method==ClpSolve::usePrimalorSprint) {
4573                // primal
4574                // if slp turn everything off
4575                if (slpValue>0) {
4576                  doCrash=false;
4577                  doSprint=0;
4578                  doIdiot=-1;
4579                  solveOptions.setSpecialOption(1,10,slpValue); // slp
4580                  method=ClpSolve::usePrimal;
4581                }
4582                if (doCrash) {
4583                  solveOptions.setSpecialOption(1,1,doCrash); // crash
4584                } else if (doSprint>0) {
4585                  // sprint overrides idiot
4586                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
4587                } else if (doIdiot>0) {
4588                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
4589                } else if (slpValue<=0) {
4590                  if (doIdiot==0) {
4591                    if (doSprint==0)
4592                      solveOptions.setSpecialOption(1,4); // all slack
4593                    else
4594                      solveOptions.setSpecialOption(1,9); // all slack or sprint
4595                  } else {
4596                    if (doSprint==0)
4597                      solveOptions.setSpecialOption(1,8); // all slack or idiot
4598                    else
4599                      solveOptions.setSpecialOption(1,7); // initiative
4600                  }
4601                }
4602                if (basisHasValues==-1)
4603                  solveOptions.setSpecialOption(1,11); // switch off values
4604              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
4605                int barrierOptions = choleskyType;
4606                if (scaleBarrier)
4607                  barrierOptions |= 8;
4608                if (doKKT)
4609                  barrierOptions |= 16;
4610                if (gamma)
4611                  barrierOptions |= 32*gamma;
4612                if (crossover==3) 
4613                  barrierOptions |= 256; // try presolve in crossover
4614                solveOptions.setSpecialOption(4,barrierOptions);
4615              }
4616              model2->setMaximumSeconds(model_.getMaximumSeconds());
4617#ifdef COIN_HAS_LINK
4618              OsiSolverInterface * coinSolver = model_.solver();
4619              OsiSolverLink * linkSolver = dynamic_cast< OsiSolverLink*> (coinSolver);
4620              if (!linkSolver) {
4621                model2->initialSolve(solveOptions);
4622              } else {
4623                // special solver
4624                int testOsiOptions = parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].intValue();
4625                double * solution = NULL;
4626                if (testOsiOptions<10) {
4627                  solution = linkSolver->nonlinearSLP(slpValue>0 ? slpValue : 20 ,1.0e-5);
4628                } else if (testOsiOptions>=10) {
4629                  CoinModel coinModel = *linkSolver->coinModel();
4630                  ClpSimplex * tempModel = approximateSolution(coinModel,slpValue>0 ? slpValue : 50 ,1.0e-5,0);
4631                  assert (tempModel);
4632                  solution = CoinCopyOfArray(tempModel->primalColumnSolution(),coinModel.numberColumns());
4633                  model2->setObjectiveValue(tempModel->objectiveValue());
4634                  model2->setProblemStatus(tempModel->problemStatus());
4635                  model2->setSecondaryStatus(tempModel->secondaryStatus());
4636                  delete tempModel;
4637                }
4638                if (solution) {
4639                  memcpy(model2->primalColumnSolution(),solution,
4640                         CoinMin(model2->numberColumns(),linkSolver->coinModel()->numberColumns())*sizeof(double));
4641                  delete [] solution;
4642                } else {
4643                  printf("No nonlinear solution\n");
4644                }
4645              }
4646#else
4647              model2->initialSolve(solveOptions);
4648#endif
4649              {
4650                // map states
4651                /* clp status
4652                   -1 - unknown e.g. before solve or if postSolve says not optimal
4653                   0 - optimal
4654                   1 - primal infeasible
4655                   2 - dual infeasible
4656                   3 - stopped on iterations or time
4657                   4 - stopped due to errors
4658                   5 - stopped by event handler (virtual int ClpEventHandler::event()) */
4659                /* cbc status
4660                   -1 before branchAndBound
4661                   0 finished - check isProvenOptimal or isProvenInfeasible to see if solution found
4662                   (or check value of best solution)
4663                   1 stopped - on maxnodes, maxsols, maxtime
4664                   2 difficulties so run was abandoned
4665                   (5 event user programmed event occurred) */
4666                /* clp secondary status of problem - may get extended
4667                   0 - none
4668                   1 - primal infeasible because dual limit reached OR probably primal
4669                   infeasible but can't prove it (main status 4)
4670                   2 - scaled problem optimal - unscaled problem has primal infeasibilities
4671                   3 - scaled problem optimal - unscaled problem has dual infeasibilities
4672                   4 - scaled problem optimal - unscaled problem has primal and dual infeasibilities
4673                   5 - giving up in primal with flagged variables
4674                   6 - failed due to empty problem check
4675                   7 - postSolve says not optimal
4676                   8 - failed due to bad element check
4677                   9 - status was 3 and stopped on time
4678                   100 up - translation of enum from ClpEventHandler
4679                */
4680                /* cbc secondary status of problem
4681                   -1 unset (status_ will also be -1)
4682                   0 search completed with solution
4683                   1 linear relaxation not feasible (or worse than cutoff)
4684                   2 stopped on gap
4685                   3 stopped on nodes
4686                   4 stopped on time
4687                   5 stopped on user event
4688                   6 stopped on solutions
4689                   7 linear relaxation unbounded
4690                */
4691                int iStatus = model2->status();
4692                int iStatus2 = model2->secondaryStatus();
4693                if (iStatus==0) {
4694                  iStatus2=0;
4695                } else if (iStatus==1) {
4696                  iStatus=0;
4697                  iStatus2=1; // say infeasible
4698                } else if (iStatus==2) {
4699                  iStatus=0;
4700                  iStatus2=7; // say unbounded
4701                } else if (iStatus==3) {
4702                  iStatus=1;
4703                  if (iStatus2==9)
4704                    iStatus2=4;
4705                  else
4706                    iStatus2=3; // Use nodes - as closer than solutions
4707                } else if (iStatus==4) {
4708                  iStatus=2; // difficulties
4709                  iStatus2=0; 
4710                }
4711                model_.setProblemStatus(iStatus);
4712                model_.setSecondaryStatus(iStatus2);
4713                assert (lpSolver==clpSolver->getModelPtr());
4714                assert (clpSolver==model_.solver());
4715                clpSolver->setWarmStart(NULL);
4716                // and in babModel if exists
4717                if (babModel_) {
4718                  babModel_->setProblemStatus(iStatus);
4719                  babModel_->setSecondaryStatus(iStatus2);
4720                } 
4721#ifdef NEW_STYLE_SOLVER
4722                int returnCode = callBack_->callBack(&model_,1);
4723#else
4724                int returnCode=callBack(&model,1);
4725#endif
4726                if (returnCode) {
4727                  // exit if user wants
4728#ifdef NEW_STYLE_SOLVER
4729                  updateModel(model2,returnMode);
4730#else
4731                  delete babModel_;
4732                  babModel_ = NULL;
4733#endif
4734                  return returnCode;
4735                }
4736              }
4737              basisHasValues=1;
4738              if (dualize) {
4739                int returnCode=((ClpSimplexOther *) lpSolver)->restoreFromDual(model2);
4740                if (model2->status()==3)
4741                  returnCode=0;
4742                delete model2;
4743                if (returnCode&&dualize!=2)
4744                  lpSolver->primal(1);
4745                model2=lpSolver;
4746              }
4747#ifdef NEW_STYLE_SOLVER
4748              updateModel(model2,returnMode);
4749              for (iUser=0;iUser<numberUserFunctions_;iUser++) {
4750                if (statusUserFunction_[iUser])
4751                  userFunction_[iUser]->exportSolution(this,1);
4752              }
4753#endif
4754#ifdef COIN_HAS_ASL
4755              if (statusUserFunction_[0]) {
4756                double value = model2->getObjValue()*model2->getObjSense();
4757                char buf[300];
4758                int pos=0;
4759                int iStat = model2->status();
4760                if (iStat==0) {
4761                  pos += sprintf(buf+pos,"optimal," );
4762                } else if (iStat==1) {
4763                  // infeasible
4764                  pos += sprintf(buf+pos,"infeasible,");
4765                } else if (iStat==2) {
4766                  // unbounded
4767                  pos += sprintf(buf+pos,"unbounded,");
4768                } else if (iStat==3) {
4769                  pos += sprintf(buf+pos,"stopped on iterations or time,");
4770                } else if (iStat==4) {
4771                  iStat = 7;
4772                  pos += sprintf(buf+pos,"stopped on difficulties,");
4773                } else if (iStat==5) {
4774                  iStat = 3;
4775                  pos += sprintf(buf+pos,"stopped on ctrl-c,");
4776                } else {
4777                  pos += sprintf(buf+pos,"status unknown,");
4778                  iStat=6;
4779                }
4780                info.problemStatus=iStat;
4781                info.objValue = value;
4782                pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
4783                               value);
4784                sprintf(buf+pos,"\n%d iterations",
4785                        model2->getIterationCount());
4786                free(info.primalSolution);
4787                int numberColumns=model2->numberColumns();
4788                info.primalSolution = (double *) malloc(numberColumns*sizeof(double));
4789                CoinCopyN(model2->primalColumnSolution(),numberColumns,info.primalSolution);
4790                int numberRows = model2->numberRows();
4791                free(info.dualSolution);
4792                info.dualSolution = (double *) malloc(numberRows*sizeof(double));
4793                CoinCopyN(model2->dualRowSolution(),numberRows,info.dualSolution);
4794                CoinWarmStartBasis * basis = model2->getBasis();
4795                free(info.rowStatus);
4796                info.rowStatus = (int *) malloc(numberRows*sizeof(int));
4797                free(info.columnStatus);
4798                info.columnStatus = (int *) malloc(numberColumns*sizeof(int));
4799                // Put basis in
4800                int i;
4801                // free,basic,ub,lb are 0,1,2,3
4802                for (i=0;i<numberRows;i++) {
4803                  CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
4804                  info.rowStatus[i]=status;
4805                }
4806                for (i=0;i<numberColumns;i++) {
4807                  CoinWarmStartBasis::Status status = basis->getStructStatus(i);
4808                  info.columnStatus[i]=status;
4809                }
4810                // put buffer into info
4811                strcpy(info.buffer,buf);
4812                delete basis;
4813              }
4814#endif
4815            } else {
4816#ifndef DISALLOW_PRINTING
4817              std::cout<<"** Current model not valid"<<std::endl;
4818#endif
4819            }
4820            break;
4821          case STATISTICS:
4822            if (goodModel) {
4823              // If presolve on look at presolved
4824              bool deleteModel2=false;
4825              ClpSimplex * model2 = lpSolver;
4826              if (preSolve) {
4827                ClpPresolve pinfo;
4828                int presolveOptions2 = presolveOptions&~0x40000000;
4829                if ((presolveOptions2&0xffff)!=0)
4830                  pinfo.setPresolveActions(presolveOptions2);
4831                pinfo.setSubstitution(substitution);
4832                if ((printOptions&1)!=0)
4833                  pinfo.statistics();
4834                double presolveTolerance = 
4835                  parameters_[whichParam(PRESOLVETOLERANCE,numberParameters_,parameters_)].doubleValue();
4836                model2 = 
4837                  pinfo.presolvedModel(*lpSolver,presolveTolerance,
4838                                       true,preSolve);
4839                if (model2) {
4840                  printf("Statistics for presolved model\n");
4841                  deleteModel2=true;
4842                } else {
4843                  printf("Presolved model looks infeasible - will use unpresolved\n");
4844                  model2 = lpSolver;
4845                }
4846              } else {
4847                printf("Statistics for unpresolved model\n");
4848                model2 =  lpSolver;
4849              }
4850              statistics(lpSolver,model2);
4851              if (deleteModel2)
4852                delete model2;
4853            } else {
4854#ifndef DISALLOW_PRINTING
4855              std::cout<<"** Current model not valid"<<std::endl;
4856#endif
4857            }
4858            break;
4859          case TIGHTEN:
4860            if (goodModel) {
4861              int numberInfeasibilities = lpSolver->tightenPrimalBounds();
4862              if (numberInfeasibilities)
4863                std::cout<<"** Analysis indicates model infeasible"<<std::endl;
4864            } else {
4865#ifndef DISALLOW_PRINTING
4866              std::cout<<"** Current model not valid"<<std::endl;
4867#endif
4868            }
4869            break;
4870          case PLUSMINUS:
4871            if (goodModel) {
4872              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
4873              ClpPackedMatrix* clpMatrix =
4874                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
4875              if (clpMatrix) {
4876                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
4877                if (newMatrix->getIndices()) {
4878                  lpSolver->replaceMatrix(newMatrix);
4879                  delete saveMatrix;
4880                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
4881                } else {
4882                  std::cout<<"Matrix can not be converted to +- 1 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 OUTDUPROWS:
4894            if (goodModel) {
4895              int numberRows = clpSolver->getNumRows();
4896              //int nOut = outDupRow(clpSolver);
4897              CglDuplicateRow dupcuts(clpSolver);
4898              storedCuts = dupcuts.outDuplicates(clpSolver)!=0;
4899              int nOut = numberRows-clpSolver->getNumRows();
4900              if (nOut&&!noPrinting_)
4901                sprintf(generalPrint,"%d rows eliminated",nOut);
4902              generalMessageHandler->message(CLP_GENERAL,generalMessages)
4903                << generalPrint
4904                <<CoinMessageEol;
4905            } else {
4906#ifndef DISALLOW_PRINTING
4907              std::cout<<"** Current model not valid"<<std::endl;
4908#endif
4909            }
4910            break;
4911          case NETWORK:
4912            if (goodModel) {
4913              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
4914              ClpPackedMatrix* clpMatrix =
4915                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
4916              if (clpMatrix) {
4917                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
4918                if (newMatrix->getIndices()) {
4919                  lpSolver->replaceMatrix(newMatrix);
4920                  delete saveMatrix;
4921                  std::cout<<"Matrix converted to network matrix"<<std::endl;
4922                } else {
4923                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
4924                }
4925              } else {
4926                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
4927              }
4928            } else {
4929#ifndef DISALLOW_PRINTING
4930              std::cout<<"** Current model not valid"<<std::endl;
4931#endif
4932            }
4933            break;
4934          case DOHEURISTIC:
4935            if (goodModel) {
4936              int vubAction = parameters_[whichParam(VUBTRY,numberParameters_,parameters_)].intValue();
4937              if (vubAction!=-1) {
4938                // look at vubs
4939                // Just ones which affect >= extra3
4940                int extra3 = parameters_[whichParam(EXTRA3,numberParameters_,parameters_)].intValue();
4941                /* 3 is fraction of integer variables fixed (0.97)
4942                   4 is fraction of all variables fixed (0.0)
4943                */
4944                double dextra[5];
4945                dextra[1] = parameters_[whichParam(DEXTRA1,numberParameters_,parameters_)].doubleValue();
4946                dextra[2] = parameters_[whichParam(DEXTRA2,numberParameters_,parameters_)].doubleValue();
4947                dextra[3] = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
4948                dextra[4] = parameters_[whichParam(DEXTRA4,numberParameters_,parameters_)].doubleValue();
4949                if (!dextra[3])
4950                  dextra[3] = 0.97;
4951                OsiClpSolverInterface * newSolver = fixVubs(model_,extra3,vubAction,generalMessageHandler,
4952                                                            debugValues,dextra);
4953                assert (!newSolver);
4954              }
4955              // Actually do heuristics
4956              doHeuristics(&model_,2);
4957              if (model_.bestSolution()) {
4958                model_.setProblemStatus(1);
4959                model_.setSecondaryStatus(6);
4960              }
4961#ifdef NEW_STYLE_SOLVER
4962              int returnCode = callBack_->callBack(&model_,6);
4963#else
4964              int returnCode=callBack(&model,6);
4965#endif
4966              if (returnCode) {
4967                // exit if user wants
4968#ifdef NEW_STYLE_SOLVER
4969                updateModel(NULL,returnMode);
4970#else
4971                delete babModel_;
4972                babModel_ = NULL;
4973#endif
4974                return returnCode;
4975              }
4976            }
4977            break;
4978          case MIPLIB:
4979            // User can set options - main difference is lack of model and CglPreProcess
4980            goodModel=true;
4981/*
4982  Run branch-and-cut. First set a few options -- node comparison, scaling.
4983  Print elapsed time at the end.
4984*/
4985          case BAB: // branchAndBound
4986          case STRENGTHEN:
4987            if (goodModel) {
4988              bool miplib = type==MIPLIB;
4989              int logLevel = parameters_[slog].intValue();
4990              // Reduce printout
4991              if (logLevel<=1)
4992                model_.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
4993              {
4994                OsiSolverInterface * solver = model_.solver();
4995                OsiClpSolverInterface * si =
4996                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
4997                assert (si != NULL);
4998                // See if quadratic
4999#ifdef COIN_HAS_LINK
5000                if (!complicatedInteger) {
5001                  ClpSimplex * lpSolver = si->getModelPtr();
5002                  ClpQuadraticObjective * obj = (dynamic_cast< ClpQuadraticObjective*>(lpSolver->objectiveAsObject()));
5003                  if (obj) {
5004                    preProcess=0;
5005                    int testOsiOptions = parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].intValue();
5006                    parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].setIntValue(CoinMax(0,testOsiOptions));
5007                    // create coin model
5008                    coinModel = lpSolver->createCoinModel();
5009                    assert (coinModel);
5010                    // load from coin model
5011                    OsiSolverLink solver1;
5012                    OsiSolverInterface * solver2 = solver1.clone();
5013                    model_.assignSolver(solver2,false);
5014                    OsiSolverLink * si =
5015                      dynamic_cast<OsiSolverLink *>(model_.solver()) ;
5016                    assert (si != NULL);
5017                    si->setDefaultMeshSize(0.001);
5018                    // need some relative granularity
5019                    si->setDefaultBound(100.0);
5020                    double dextra3 = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
5021                    if (dextra3)
5022                      si->setDefaultMeshSize(dextra3);
5023                    si->setDefaultBound(1000.0);
5024                    si->setIntegerPriority(1000);
5025                    si->setBiLinearPriority(10000);
5026                    si->setSpecialOptions2(2+4+8);
5027                    CoinModel * model2 = (CoinModel *) coinModel;
5028                    si->load(*model2,true, parameters_[log].intValue());
5029                    // redo
5030                    solver = model_.solver();
5031                    clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
5032                    lpSolver = clpSolver->getModelPtr();
5033                    clpSolver->messageHandler()->setLogLevel(0) ;
5034                    testOsiParameters=0;
5035                    complicatedInteger=2;  // allow cuts
5036                    OsiSolverInterface * coinSolver = model_.solver();
5037                    OsiSolverLink * linkSolver = dynamic_cast< OsiSolverLink*> (coinSolver);
5038                    if (linkSolver->quadraticModel()) {
5039                      ClpSimplex * qp = linkSolver->quadraticModel();
5040                      //linkSolver->nonlinearSLP(CoinMax(slpValue,10),1.0e-5);
5041                      qp->nonlinearSLP(CoinMax(slpValue,40),1.0e-5);
5042                      qp->primal(1);
5043                      OsiSolverLinearizedQuadratic solver2(qp);
5044                      const double * solution=NULL;
5045                      // Reduce printout
5046                      solver2.setHintParam(OsiDoReducePrint,true,OsiHintTry);
5047                      CbcModel model2(solver2);
5048                      // Now do requested saves and modifications
5049                      CbcModel * cbcModel = & model2;
5050                      OsiSolverInterface * osiModel = model2.solver();
5051                      OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
5052                      ClpSimplex * clpModel = osiclpModel->getModelPtr();
5053                     
5054                      // Set changed values
5055                     
5056                      CglProbing probing;
5057                      probing.setMaxProbe(10);
5058                      probing.setMaxLook(10);
5059                      probing.setMaxElements(200);
5060                      probing.setMaxProbeRoot(50);
5061                      probing.setMaxLookRoot(10);
5062                      probing.setRowCuts(3);
5063                      probing.setUsingObjective(true);
5064                      cbcModel->addCutGenerator(&probing,-1,"Probing",true,false,false,-100,-1,-1);
5065                      cbcModel->cutGenerator(0)->setTiming(true);
5066                     
5067                      CglGomory gomory;
5068                      gomory.setLimitAtRoot(512);
5069                      cbcModel->addCutGenerator(&gomory,-98,"Gomory",true,false,false,-100,-1,-1);
5070                      cbcModel->cutGenerator(1)->setTiming(true);
5071                     
5072                      CglKnapsackCover knapsackCover;
5073                      cbcModel->addCutGenerator(&knapsackCover,-98,"KnapsackCover",true,false,false,-100,-1,-1);
5074                      cbcModel->cutGenerator(2)->setTiming(true);
5075                     
5076                      CglRedSplit redSplit;
5077                      cbcModel->addCutGenerator(&redSplit,-99,"RedSplit",true,false,false,-100,-1,-1);
5078                      cbcModel->cutGenerator(3)->setTiming(true);
5079                     
5080                      CglClique clique;
5081                      clique.setStarCliqueReport(false);
5082                      clique.setRowCliqueReport(false);
5083                      clique.setMinViolation(0.1);
5084                      cbcModel->addCutGenerator(&clique,-98,"Clique",true,false,false,-100,-1,-1);
5085                      cbcModel->cutGenerator(4)->setTiming(true);
5086                     
5087                      CglMixedIntegerRounding2 mixedIntegerRounding2;
5088                      cbcModel->addCutGenerator(&mixedIntegerRounding2,-98,"MixedIntegerRounding2",true,false,false,-100,-1,-1);
5089                      cbcModel->cutGenerator(5)->setTiming(true);
5090                     
5091                      CglFlowCover flowCover;
5092                      cbcModel->addCutGenerator(&flowCover,-98,"FlowCover",true,false,false,-100,-1,-1);
5093                      cbcModel->cutGenerator(6)->setTiming(true);
5094                     
5095                      CglTwomir twomir;
5096                      twomir.setMaxElements(250);
5097                      cbcModel->addCutGenerator(&twomir,-99,"Twomir",true,false,false,-100,-1,-1);
5098                      cbcModel->cutGenerator(7)->setTiming(true);
5099                     
5100                      CbcHeuristicFPump heuristicFPump(*cbcModel);
5101                      heuristicFPump.setWhen(13);
5102                      heuristicFPump.setMaximumPasses(20);
5103                      heuristicFPump.setMaximumRetries(7);
5104                      heuristicFPump.setHeuristicName("feasibility pump");
5105                      heuristicFPump.setInitialWeight(1);
5106                      heuristicFPump.setFractionSmall(0.6);
5107                      cbcModel->addHeuristic(&heuristicFPump);
5108                     
5109                      CbcRounding rounding(*cbcModel);
5110                      rounding.setHeuristicName("rounding");
5111                      cbcModel->addHeuristic(&rounding);
5112                     
5113                      CbcHeuristicLocal heuristicLocal(*cbcModel);
5114                      heuristicLocal.setHeuristicName("combine solutions");
5115                      heuristicLocal.setSearchType(1);
5116                      heuristicLocal.setFractionSmall(0.6);
5117                      cbcModel->addHeuristic(&heuristicLocal);
5118                     
5119                      CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
5120                      heuristicGreedyCover.setHeuristicName("greedy cover");
5121                      cbcModel->addHeuristic(&heuristicGreedyCover);
5122                     
5123                      CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
5124                      heuristicGreedyEquality.setHeuristicName("greedy equality");
5125                      cbcModel->addHeuristic(&heuristicGreedyEquality);
5126                     
5127                      CbcCompareDefault compare;
5128                      cbcModel->setNodeComparison(compare);
5129                      cbcModel->setNumberBeforeTrust(5);
5130                      cbcModel->setSpecialOptions(2);
5131                      cbcModel->messageHandler()->setLogLevel(1);
5132                      cbcModel->setMaximumCutPassesAtRoot(-100);
5133                      cbcModel->setMaximumCutPasses(1);
5134                      cbcModel->setMinimumDrop(0.05);
5135                      // For branchAndBound this may help
5136                      clpModel->defaultFactorizationFrequency();
5137                      clpModel->setDualBound(1.0001e+08);
5138                      clpModel->setPerturbation(50);
5139                      osiclpModel->setSpecialOptions(193);
5140                      osiclpModel->messageHandler()->setLogLevel(0);
5141                      osiclpModel->setIntParam(OsiMaxNumIterationHotStart,100);
5142                      osiclpModel->setHintParam(OsiDoReducePrint,true,OsiHintTry);
5143                      // You can save some time by switching off message building
5144                      // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
5145                     
5146                      // Solve
5147                     
5148                      cbcModel->initialSolve();
5149                      if (clpModel->tightenPrimalBounds()!=0) {
5150#ifndef DISALLOW_PRINTING
5151                        std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
5152#endif
5153                        break;
5154                      }
5155                      clpModel->dual();  // clean up
5156                      cbcModel->initialSolve();
5157#ifdef CBC_THREAD
5158                      int numberThreads =parameters_[whichParam(THREADS,numberParameters_,parameters_)].intValue();
5159                      cbcModel->setNumberThreads(numberThreads%100);
5160                      cbcModel->setThreadMode(numberThreads/100);
5161#endif
5162                      cbcModel->branchAndBound();
5163                      OsiSolverLinearizedQuadratic * solver3 = dynamic_cast<OsiSolverLinearizedQuadratic *> (model2.solver());
5164                      assert (solver3);
5165                      solution = solver3->bestSolution();
5166                      double bestObjectiveValue = solver3->bestObjectiveValue();
5167                      linkSolver->setBestObjectiveValue(bestObjectiveValue);
5168                      linkSolver->setBestSolution(solution,solver3->getNumCols());
5169                      CbcHeuristicDynamic3 dynamic(model_);
5170                      dynamic.setHeuristicName("dynamic pass thru");
5171                      model_.addHeuristic(&dynamic);
5172                      // if convex
5173                      if ((linkSolver->specialOptions2()&4)!=0) {
5174                        int numberColumns = coinModel->numberColumns();
5175                        assert (linkSolver->objectiveVariable()==numberColumns);
5176                        // add OA cut
5177                        double offset;
5178                        double * gradient = new double [numberColumns+1];
5179                        memcpy(gradient,qp->objectiveAsObject()->gradient(qp,solution,offset,true,2),
5180                               numberColumns*sizeof(double));
5181                        double rhs = 0.0;
5182                        int * column = new int[numberColumns+1];
5183                        int n=0;
5184                        for (int i=0;i<numberColumns;i++) {
5185                          double value = gradient[i];
5186                          if (fabs(value)>1.0e-12) {
5187                            gradient[n]=value;
5188                            rhs += value*solution[i];
5189                            column[n++]=i;
5190                          }
5191                        }
5192                        gradient[n]=-1.0;
5193                        column[n++]=numberColumns;
5194                        storedAmpl.addCut(-COIN_DBL_MAX,offset+1.0e-7,n,column,gradient);
5195                        delete [] gradient;
5196                        delete [] column;
5197                      }
5198                      // could do three way branching round a) continuous b) best solution
5199                      printf("obj %g\n",bestObjectiveValue);
5200                      linkSolver->initialSolve();
5201                    }
5202                  }
5203                }
5204#endif
5205                si->setSpecialOptions(0x40000000);
5206              }
5207              if (!miplib) {
5208                if (!preSolve) {
5209                  model_.solver()->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry);
5210                  model_.solver()->setHintParam(OsiDoPresolveInResolve,false,OsiHintTry);
5211                }
5212                double time1a = CoinCpuTime();
5213                model_.initialSolve();
5214                OsiSolverInterface * solver = model_.solver();
5215                OsiClpSolverInterface * si =
5216                  dynamic_cast<OsiClpSolverInterface *>(solver) ;
5217                ClpSimplex * clpSolver = si->getModelPtr();
5218                int iStatus = clpSolver->status();
5219                int iStatus2 = clpSolver->secondaryStatus();
5220                if (iStatus==0) {
5221                  iStatus2=0;
5222                } else if (iStatus==1) {
5223                  iStatus=0;
5224                  iStatus2=1; // say infeasible
5225                } else if (iStatus==2) {
5226                  iStatus=0;
5227                  iStatus2=7; // say unbounded
5228                } else if (iStatus==3) {
5229                  iStatus=1;
5230                  if (iStatus2==9)
5231                    iStatus2=4;
5232                  else
5233                    iStatus2=3; // Use nodes - as closer than solutions
5234                } else if (iStatus==4) {
5235                  iStatus=2; // difficulties
5236                  iStatus2=0; 
5237                }
5238                model_.setProblemStatus(iStatus);
5239                model_.setSecondaryStatus(iStatus2);
5240                si->setWarmStart(NULL);
5241#ifdef NEW_STYLE_SOLVER
5242                int returnCode = callBack_->callBack(&model_,1);
5243#else
5244                int returnCode=callBack(&model_,1);
5245#endif
5246                if (returnCode) {
5247                  // exit if user wants
5248#ifdef NEW_STYLE_SOLVER
5249                  updateModel(NULL,returnMode);
5250#else
5251                  delete babModel_;
5252                  babModel_ = NULL;
5253#endif
5254                  return returnCode;
5255                }
5256                clpSolver->setSpecialOptions(clpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
5257                if (!noPrinting_) {
5258                  sprintf(generalPrint,"Continuous objective value is %g - %.2f seconds",
5259                          solver->getObjValue(),CoinCpuTime()-time1a);