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

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

correct reported value of obj when maxim

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