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

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

changes for heuristic clone

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