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

Last change on this file since 905 was 905, checked in by ladanyi, 11 years ago

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

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