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

Last change on this file since 1131 was 1131, checked in by forrest, 10 years ago

add release version

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