source: stable/2.2/Cbc/src/CbcSolver.cpp @ 1048

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

2.1 changes to 2.2 for Stefan

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