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

Last change on this file since 1020 was 1020, checked in by forrest, 11 years ago

tweaks to heuristics

File size: 353.5 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  int statusUserFunction_[1];
3624  int numberUserFunctions_=1; // to allow for ampl
3625#else
3626  delete babModel_;
3627  babModel_ = NULL;
3628  CbcOrClpRead_mode=1;
3629  delete [] statusUserFunction_;
3630  statusUserFunction_ = new int [numberUserFunctions_];
3631  int iUser;
3632#endif
3633  // Statistics
3634  double statistics_seconds=0.0, statistics_obj=0.0;
3635  double statistics_continuous=0.0, statistics_tighter=0.0;
3636  double statistics_cut_time=0.0;
3637  int statistics_nodes=0, statistics_iterations=0;
3638  int statistics_nrows=0, statistics_ncols=0;
3639  int statistics_nprocessedrows=0, statistics_nprocessedcols=0;
3640  std::string statistics_result;
3641  int * statistics_number_cuts=NULL;
3642  const char ** statistics_name_generators=NULL;
3643  int statistics_number_generators=0;
3644  memset(statusUserFunction_,0,numberUserFunctions_*sizeof(int));
3645  /* Note
3646     This is meant as a stand-alone executable to do as much of coin as possible.
3647     It should only have one solver known to it.
3648  */
3649  OsiClpSolverInterface * originalSolver = dynamic_cast<OsiClpSolverInterface *> (model_.solver());
3650  assert (originalSolver);
3651  CoinMessageHandler * generalMessageHandler = model_.messageHandler();
3652  generalMessageHandler->setPrefix(false);
3653  // Move handler across if not default
3654  if (!originalSolver->defaultHandler()&&originalSolver->getModelPtr()->defaultHandler())
3655    originalSolver->getModelPtr()->passInMessageHandler(originalSolver->messageHandler());
3656  CoinMessages generalMessages = originalSolver->getModelPtr()->messages();
3657  char generalPrint[10000];
3658  if (originalSolver->getModelPtr()->logLevel()==0)
3659    noPrinting=true;
3660#if NEW_STYLE_SOLVER==0
3661  bool noPrinting_=noPrinting;
3662#endif
3663  // see if log in list
3664  for (int i=1;i<argc;i++) {
3665    if (!strncmp(argv[i],"log",3)) {
3666      const char * equals = strchr(argv[i],'=');
3667      if (equals&&atoi(equals+1)!=0) 
3668        noPrinting_=false;
3669      else
3670        noPrinting_=true;
3671      break;
3672    } else if (!strncmp(argv[i],"-log",4)&&i<argc-1) {
3673      if (atoi(argv[i+1])!=0) 
3674        noPrinting_=false;
3675      else
3676        noPrinting_=true;
3677      break;
3678    }
3679  }
3680  double time0;
3681  {
3682    double time1 = CoinCpuTime(),time2;
3683    time0=time1;
3684    bool goodModel=(originalSolver->getNumCols()) ? true : false;
3685
3686    CoinSighandler_t saveSignal=SIG_DFL;
3687    // register signal handler
3688    saveSignal = signal(SIGINT,signal_handler);
3689    // Set up all non-standard stuff
3690    int cutPass=-1234567;
3691    int cutPassInTree=-1234567;
3692    int tunePreProcess=0;
3693    int testOsiParameters=-1;
3694    // 0 normal, 1 from ampl or MIQP etc (2 allows cuts)
3695    int complicatedInteger=0;
3696    OsiSolverInterface * solver = model_.solver();
3697    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3698    ClpSimplex * lpSolver = clpSolver->getModelPtr();
3699    if (noPrinting_) {
3700      setCbcOrClpPrinting(false);
3701      lpSolver->setLogLevel(0);
3702    }
3703    // For priorities etc
3704    int * priorities=NULL;
3705    int * branchDirection=NULL;
3706    double * pseudoDown=NULL;
3707    double * pseudoUp=NULL;
3708    double * solutionIn = NULL;
3709    int * prioritiesIn = NULL;
3710    int numberSOS = 0;
3711    int * sosStart = NULL;
3712    int * sosIndices = NULL;
3713    char * sosType = NULL;
3714    double * sosReference = NULL;
3715    int * cut=NULL;
3716    int * sosPriority=NULL;
3717    CglStored storedAmpl;
3718    CoinModel * coinModel = NULL;
3719    CoinModel saveCoinModel;
3720    CoinModel saveTightenedModel;
3721    int * whichColumn = NULL;
3722    int * knapsackStart=NULL;
3723    int * knapsackRow=NULL;
3724    int numberKnapsack=0;
3725#if NEW_STYLE_SOLVER
3726    int numberInputs=0;
3727    readMode_=CbcOrClpRead_mode;
3728    for (iUser=0;iUser<numberUserFunctions_;iUser++) {
3729      int status = userFunction_[iUser]->importData(this,argc,const_cast<char **>(argv));
3730      if (status>=0) {
3731        if (!status) {
3732          numberInputs++;
3733          statusUserFunction_[iUser]=1;
3734          goodModel=true;
3735          solver = model_.solver();
3736          clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3737          lpSolver = clpSolver->getModelPtr();
3738        } else {
3739          printf("Bad input from user function %s\n",userFunction_[iUser]->name().c_str());
3740          abort();
3741        }
3742      }
3743    }
3744    if (numberInputs>1) {
3745      printf("Two or more user inputs!\n");
3746      abort();
3747    }
3748    if (originalCoinModel_)
3749      complicatedInteger=1;
3750    testOsiParameters = intValue(TESTOSI);
3751    if (noPrinting_) {
3752      model_.messageHandler()->setLogLevel(0);
3753      setCbcOrClpPrinting(false);
3754    }
3755    CbcOrClpRead_mode=readMode_;
3756#endif
3757#ifdef COIN_HAS_ASL
3758    ampl_info info;
3759    {
3760    memset(&info,0,sizeof(info));
3761    if (argc>2&&!strcmp(argv[2],"-AMPL")) {
3762      statusUserFunction_[0]=1;
3763      // see if log in list
3764      noPrinting_=true;
3765      for (int i=1;i<argc;i++) {
3766        if (!strncmp(argv[i],"log",3)) {
3767          const char * equals = strchr(argv[i],'=');
3768          if (equals&&atoi(equals+1)>0) {
3769            noPrinting_=false;
3770            info.logLevel=atoi(equals+1);
3771            int log = whichParam(LOGLEVEL,numberParameters_,parameters_);
3772            parameters_[log].setIntValue(info.logLevel);
3773            // mark so won't be overWritten
3774            info.numberRows=-1234567;
3775            break;
3776          }
3777        }
3778      }
3779
3780      union { void * voidModel; CoinModel * model; } coinModelStart;
3781      coinModelStart.model=NULL;
3782      int returnCode = readAmpl(&info,argc, const_cast<char **>(argv),& coinModelStart.voidModel);
3783      coinModel=coinModelStart.model;
3784      if (returnCode)
3785        return returnCode;
3786      CbcOrClpRead_mode=2; // so will start with parameters
3787      // see if log in list (including environment)
3788      for (int i=1;i<info.numberArguments;i++) {
3789        if (!strcmp(info.arguments[i],"log")) {
3790          if (i<info.numberArguments-1&&atoi(info.arguments[i+1])>0)
3791            noPrinting_=false;
3792          break;
3793        }
3794      }
3795      if (noPrinting_) {
3796        model_.messageHandler()->setLogLevel(0);
3797        setCbcOrClpPrinting(false);
3798      }
3799      if (!noPrinting_)
3800        printf("%d rows, %d columns and %d elements\n",
3801               info.numberRows,info.numberColumns,info.numberElements);
3802#ifdef COIN_HAS_LINK
3803      if (!coinModel) {
3804#endif
3805      solver->loadProblem(info.numberColumns,info.numberRows,info.starts,
3806                          info.rows,info.elements,
3807                          info.columnLower,info.columnUpper,info.objective,
3808                          info.rowLower,info.rowUpper);
3809      // take off cuts if ampl wants that
3810      if (info.cut&&0) {
3811        printf("AMPL CUTS OFF until global cuts fixed\n");
3812        info.cut=NULL;
3813      }
3814      if (info.cut) {
3815        int numberRows = info.numberRows;
3816        int * whichRow = new int [numberRows];
3817        // Row copy
3818        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
3819        const double * elementByRow = matrixByRow->getElements();
3820        const int * column = matrixByRow->getIndices();
3821        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
3822        const int * rowLength = matrixByRow->getVectorLengths();
3823       
3824        const double * rowLower = solver->getRowLower();
3825        const double * rowUpper = solver->getRowUpper();
3826        int nDelete=0;
3827        for (int iRow=0;iRow<numberRows;iRow++) {
3828          if (info.cut[iRow]) {
3829            whichRow[nDelete++]=iRow;
3830            int start = rowStart[iRow];
3831            storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
3832                          rowLength[iRow],column+start,elementByRow+start);
3833          }
3834        }
3835        solver->deleteRows(nDelete,whichRow);
3836        delete [] whichRow;
3837      }
3838#ifdef COIN_HAS_LINK
3839      } else {
3840        // save
3841        saveCoinModel = *coinModel;
3842        // load from coin model
3843        OsiSolverLink solver1;
3844        OsiSolverInterface * solver2 = solver1.clone();
3845        model_.assignSolver(solver2,false);
3846        OsiSolverLink * si =
3847          dynamic_cast<OsiSolverLink *>(model_.solver()) ;
3848        assert (si != NULL);
3849        si->setDefaultMeshSize(0.001);
3850        // need some relative granularity
3851        si->setDefaultBound(100.0);
3852        double dextra3 = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
3853        if (dextra3)
3854          si->setDefaultMeshSize(dextra3);
3855        si->setDefaultBound(100000.0);
3856        si->setIntegerPriority(1000);
3857        si->setBiLinearPriority(10000);
3858        CoinModel * model2 = (CoinModel *) coinModel;
3859        int logLevel = parameters_[whichParam(LOGLEVEL,numberParameters_,parameters_)].intValue();
3860        si->load(*model2,true,logLevel);
3861        // redo
3862        solver = model_.solver();
3863        clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3864        lpSolver = clpSolver->getModelPtr();
3865        clpSolver->messageHandler()->setLogLevel(0) ;
3866        testOsiParameters=0;
3867        parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].setIntValue(0);
3868        complicatedInteger=1;
3869        if (info.cut) {
3870          printf("Sorry - can't do cuts with LOS as ruins delicate row order\n");
3871          abort();
3872          int numberRows = info.numberRows;
3873          int * whichRow = new int [numberRows];
3874          // Row copy
3875          const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
3876          const double * elementByRow = matrixByRow->getElements();
3877          const int * column = matrixByRow->getIndices();
3878          const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
3879          const int * rowLength = matrixByRow->getVectorLengths();
3880         
3881          const double * rowLower = solver->getRowLower();
3882          const double * rowUpper = solver->getRowUpper();
3883          int nDelete=0;
3884          for (int iRow=0;iRow<numberRows;iRow++) {
3885            if (info.cut[iRow]) {
3886              whichRow[nDelete++]=iRow;
3887              int start = rowStart[iRow];
3888              storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
3889                                rowLength[iRow],column+start,elementByRow+start);
3890            }
3891          }
3892          solver->deleteRows(nDelete,whichRow);
3893          // and special matrix
3894          si->cleanMatrix()->deleteRows(nDelete,whichRow);
3895          delete [] whichRow;
3896        }
3897      }
3898#endif
3899      // If we had a solution use it
3900      if (info.primalSolution) {
3901        solver->setColSolution(info.primalSolution);
3902      }
3903      // status
3904      if (info.rowStatus) {
3905        unsigned char * statusArray = lpSolver->statusArray();
3906        int i;
3907        for (i=0;i<info.numberColumns;i++)
3908          statusArray[i]=(char)info.columnStatus[i];
3909        statusArray+=info.numberColumns;
3910        for (i=0;i<info.numberRows;i++)
3911          statusArray[i]=(char)info.rowStatus[i];
3912        CoinWarmStartBasis * basis = lpSolver->getBasis();
3913        solver->setWarmStart(basis);
3914        delete basis;
3915      }
3916      freeArrays1(&info);
3917      // modify objective if necessary
3918      solver->setObjSense(info.direction);
3919      solver->setDblParam(OsiObjOffset,info.offset);
3920      if (info.offset) {
3921        sprintf(generalPrint,"Ampl objective offset is %g",
3922                info.offset);
3923        generalMessageHandler->message(CLP_GENERAL,generalMessages)
3924          << generalPrint
3925          <<CoinMessageEol;
3926      }
3927      // Set integer variables (unless nonlinear when set)
3928      if (!info.nonLinear) {
3929        for (int i=info.numberColumns-info.numberIntegers;
3930             i<info.numberColumns;i++)
3931          solver->setInteger(i);
3932      }
3933      goodModel=true;
3934      // change argc etc
3935      argc = info.numberArguments;
3936      argv = const_cast<const char **>(info.arguments);
3937    }
3938    }
3939#endif   
3940    // default action on import
3941    int allowImportErrors=0;
3942    int keepImportNames=1;
3943    int doIdiot=-1;
3944    int outputFormat=2;
3945    int slpValue=-1;
3946    int cppValue=-1;
3947    int printOptions=0;
3948    int printMode=0;
3949    int presolveOptions=0;
3950    int substitution=3;
3951    int dualize=3;
3952    int doCrash=0;
3953    int doVector=0;
3954    int doSprint=-1;
3955    int doScaling=4;
3956    // set reasonable defaults
3957    int preSolve=5;
3958    int preProcess=4;
3959    bool useStrategy=false;
3960    bool preSolveFile=false;
3961    bool strongChanged=false;
3962   
3963    double djFix=1.0e100;
3964    double gapRatio=1.0e100;
3965    double tightenFactor=0.0;
3966    const char dirsep =  CoinFindDirSeparator();
3967    std::string directory;
3968    std::string dirSample;
3969    std::string dirNetlib;
3970    std::string dirMiplib;
3971    if (dirsep == '/') {
3972      directory = "./";
3973      dirSample = "../../Data/Sample/";
3974      dirNetlib = "../../Data/Netlib/";
3975      dirMiplib = "../../Data/miplib3/";
3976    } else {
3977      directory = ".\\";
3978      dirSample = "..\\..\\Data\\Sample\\";
3979      dirNetlib = "..\\..\\Data\\Netlib\\";
3980      dirMiplib = "..\\..\\Data\\miplib3\\";
3981    }
3982    std::string defaultDirectory = directory;
3983    std::string importFile ="";
3984    std::string exportFile ="default.mps";
3985    std::string importBasisFile ="";
3986    std::string importPriorityFile ="";
3987    std::string debugFile="";
3988    std::string printMask="";
3989    double * debugValues = NULL;
3990    int numberDebugValues = -1;
3991    int basisHasValues=0;
3992    std::string exportBasisFile ="default.bas";
3993    std::string saveFile ="default.prob";
3994    std::string restoreFile ="default.prob";
3995    std::string solutionFile ="stdout";
3996    std::string solutionSaveFile ="solution.file";
3997    int slog = whichParam(SOLVERLOGLEVEL,numberParameters_,parameters_);
3998    int log = whichParam(LOGLEVEL,numberParameters_,parameters_);
3999    double normalIncrement=model_.getCutoffIncrement();;
4000    if (testOsiParameters>=0) {
4001      // trying nonlinear - switch off some stuff
4002      preProcess=0;
4003    }
4004    // Set up likely cut generators and defaults
4005    int nodeStrategy=0;
4006    int doSOS=1;
4007    int verbose=0;
4008    CglGomory gomoryGen;
4009    // try larger limit
4010    gomoryGen.setLimitAtRoot(1000);
4011    gomoryGen.setLimit(50);
4012    // set default action (0=off,1=on,2=root)
4013    int gomoryAction=3;
4014
4015    CglProbing probingGen;
4016    probingGen.setUsingObjective(1);
4017    probingGen.setMaxPass(1);
4018    probingGen.setMaxPassRoot(1);
4019    // Number of unsatisfied variables to look at
4020    probingGen.setMaxProbe(10);
4021    probingGen.setMaxProbeRoot(50);
4022    // How far to follow the consequences
4023    probingGen.setMaxLook(10);
4024    probingGen.setMaxLookRoot(50);
4025    probingGen.setMaxLookRoot(10);
4026    // Only look at rows with fewer than this number of elements
4027    probingGen.setMaxElements(200);
4028    probingGen.setMaxElementsRoot(300);
4029    probingGen.setRowCuts(3);
4030    // set default action (0=off,1=on,2=root)
4031    int probingAction=1;
4032
4033    CglKnapsackCover knapsackGen;
4034    //knapsackGen.switchOnExpensive();
4035    // set default action (0=off,1=on,2=root)
4036    int knapsackAction=3;
4037
4038    CglRedSplit redsplitGen;
4039    //redsplitGen.setLimit(100);
4040    // set default action (0=off,1=on,2=root)
4041    // Off as seems to give some bad cuts
4042    int redsplitAction=0;
4043
4044    CglFakeClique cliqueGen(NULL,false);
4045    //CglClique cliqueGen(false,true);
4046    cliqueGen.setStarCliqueReport(false);
4047    cliqueGen.setRowCliqueReport(false);
4048    cliqueGen.setMinViolation(0.1);
4049    // set default action (0=off,1=on,2=root)
4050    int cliqueAction=3;
4051
4052    CglMixedIntegerRounding2 mixedGen;
4053    // set default action (0=off,1=on,2=root)
4054    int mixedAction=3;
4055
4056    CglFlowCover flowGen;
4057    // set default action (0=off,1=on,2=root)
4058    int flowAction=3;
4059
4060    CglTwomir twomirGen;
4061    twomirGen.setMaxElements(250);
4062    // set default action (0=off,1=on,2=root)
4063    int twomirAction=2;
4064#ifndef DEBUG_MALLOC
4065    CglLandP landpGen;
4066#endif
4067    // set default action (0=off,1=on,2=root)
4068    int landpAction=0;
4069    CglResidualCapacity residualCapacityGen;
4070    // set default action (0=off,1=on,2=root)
4071    int residualCapacityAction=0;
4072    // Stored cuts
4073    bool storedCuts = false;
4074
4075    int useCosts=0;
4076    // don't use input solution
4077    int useSolution=-1;
4078   
4079    // total number of commands read
4080    int numberGoodCommands=0;
4081    // Set false if user does anything advanced
4082    bool defaultSettings=true;
4083
4084    // Hidden stuff for barrier
4085    int choleskyType = 0;
4086    int gamma=0;
4087    int scaleBarrier=0;
4088    int doKKT=0;
4089    int crossover=2; // do crossover unless quadratic
4090    // For names
4091    int lengthName = 0;
4092    std::vector<std::string> rowNames;
4093    std::vector<std::string> columnNames;
4094   
4095    std::string field;
4096    if (!noPrinting_) {
4097      sprintf(generalPrint,"Coin Cbc and Clp Solver version %s, build %s",
4098              CBCVERSION,__DATE__);
4099      generalMessageHandler->message(CLP_GENERAL,generalMessages)
4100        << generalPrint
4101        <<CoinMessageEol;
4102      // Print command line
4103      if (argc>1) {
4104        sprintf(generalPrint,"command line - ");
4105        for (int i=0;i<argc;i++) {
4106          if (!argv[i])
4107            break;
4108          sprintf(generalPrint+strlen(generalPrint),"%s ",argv[i]);
4109        }
4110        generalMessageHandler->message(CLP_GENERAL,generalMessages)
4111          << generalPrint
4112          <<CoinMessageEol;
4113      }
4114    }
4115    while (1) {
4116      // next command
4117      field=CoinReadGetCommand(argc,argv);
4118      // Reset time
4119      time1 = CoinCpuTime();
4120      // adjust field if has odd trailing characters
4121      char temp [200];
4122      strcpy(temp,field.c_str());
4123      int length = strlen(temp);
4124      for (int k=length-1;k>=0;k--) {
4125        if (temp[k]<' ')
4126          length--;
4127        else
4128          break;
4129      }
4130      temp[length]='\0';
4131      field=temp;
4132      // exit if null or similar
4133      if (!field.length()) {
4134        if (numberGoodCommands==1&&goodModel) {
4135          // we just had file name - do branch and bound
4136          field="branch";
4137        } else if (!numberGoodCommands) {
4138          // let's give the sucker a hint
4139          std::cout
4140            <<"CoinSolver takes input from arguments ( - switches to stdin)"
4141            <<std::endl
4142            <<"Enter ? for list of commands or help"<<std::endl;
4143          field="-";
4144        } else {
4145          break;
4146        }
4147      }
4148     
4149      // see if ? at end
4150      int numberQuery=0;
4151      if (field!="?"&&field!="???") {
4152        int length = field.length();
4153        int i;
4154        for (i=length-1;i>0;i--) {
4155          if (field[i]=='?') 
4156            numberQuery++;
4157          else
4158            break;
4159        }
4160        field=field.substr(0,length-numberQuery);
4161      }
4162      // find out if valid command
4163      int iParam;
4164      int numberMatches=0;
4165      int firstMatch=-1;
4166      for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4167        int match = parameters_[iParam].matches(field);
4168        if (match==1) {
4169          numberMatches = 1;
4170          firstMatch=iParam;
4171          break;
4172        } else {
4173          if (match&&firstMatch<0)
4174            firstMatch=iParam;
4175          numberMatches += match>>1;
4176        }
4177      }
4178      if (iParam<numberParameters_&&!numberQuery) {
4179        // found
4180        CbcOrClpParam found = parameters_[iParam];
4181        CbcOrClpParameterType type = found.type();
4182        int valid;
4183        numberGoodCommands++;
4184        if (type==BAB&&goodModel) {
4185          // check if any integers
4186#ifdef COIN_HAS_ASL
4187          if (info.numberSos&&doSOS&&statusUserFunction_[0]) {
4188            // SOS
4189            numberSOS = info.numberSos;
4190          }
4191#endif
4192          lpSolver = clpSolver->getModelPtr();
4193          if (!lpSolver->integerInformation()&&!numberSOS&&
4194              !clpSolver->numberSOS()&&!model_.numberObjects()&&!clpSolver->numberObjects())
4195            type=DUALSIMPLEX;
4196        }
4197        if (type==GENERALQUERY) {
4198          bool evenHidden=false;
4199          if ((verbose&8)!=0) {
4200            // even hidden
4201            evenHidden = true;
4202            verbose &= ~8;
4203          }
4204#ifdef COIN_HAS_ASL
4205          if (verbose<4&&statusUserFunction_[0])
4206            verbose +=4;
4207#endif
4208          if (verbose<4) {
4209            std::cout<<"In argument list keywords have leading - "
4210              ", -stdin or just - switches to stdin"<<std::endl;
4211            std::cout<<"One command per line (and no -)"<<std::endl;
4212            std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
4213            std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
4214            std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
4215            std::cout<<"abcd value sets value"<<std::endl;
4216            std::cout<<"Commands are:"<<std::endl;
4217          } else {
4218            std::cout<<"Cbc options are set within AMPL with commands like:"<<std::endl<<std::endl;
4219            std::cout<<"         option cbc_options \"cuts=root log=2 feas=on slog=1\""<<std::endl<<std::endl;
4220            std::cout<<"only maximize, dual, primal, help and quit are recognized without ="<<std::endl;
4221          }
4222          int maxAcross=5;
4223          if ((verbose%4)!=0)
4224            maxAcross=1;
4225          int limits[]={1,51,101,151,201,251,301,351,401};
4226          std::vector<std::string> types;
4227          types.push_back("Double parameters:");
4228          types.push_back("Branch and Cut double parameters:");
4229          types.push_back("Integer parameters:");
4230          types.push_back("Branch and Cut integer parameters:");
4231          types.push_back("Keyword parameters:");
4232          types.push_back("Branch and Cut keyword parameters:");
4233          types.push_back("Actions or string parameters:");
4234          types.push_back("Branch and Cut actions:");
4235          int iType;
4236          for (iType=0;iType<8;iType++) {
4237            int across=0;
4238            if ((verbose%4)!=0)
4239              std::cout<<std::endl;
4240            std::cout<<types[iType]<<std::endl;
4241            if ((verbose&2)!=0)
4242              std::cout<<std::endl;
4243            for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4244              int type = parameters_[iParam].type();
4245              if ((parameters_[iParam].displayThis()||evenHidden)&&
4246                  type>=limits[iType]
4247                  &&type<limits[iType+1]) {
4248                // but skip if not useful for ampl (and in ampl mode)
4249                if (verbose>=4&&(parameters_[iParam].whereUsed()&4)==0)
4250                  continue;
4251                if (!across) {
4252                  if ((verbose&2)==0) 
4253                    std::cout<<"  ";
4254                  else
4255                    std::cout<<"Command ";
4256                }
4257                std::cout<<parameters_[iParam].matchName()<<"  ";
4258                across++;
4259                if (across==maxAcross) {
4260                  across=0;
4261                  if ((verbose%4)!=0) {
4262                    // put out description as well
4263                    if ((verbose&1)!=0) 
4264                      std::cout<<parameters_[iParam].shortHelp();
4265                    std::cout<<std::endl;
4266                    if ((verbose&2)!=0) {
4267                      std::cout<<"---- description"<<std::endl;
4268                      parameters_[iParam].printLongHelp();
4269                      std::cout<<"----"<<std::endl<<std::endl;
4270                    }
4271                  } else {
4272                    std::cout<<std::endl;
4273                  }
4274                }
4275              }
4276            }
4277            if (across)
4278              std::cout<<std::endl;
4279          }
4280        } else if (type==FULLGENERALQUERY) {
4281          std::cout<<"Full list of commands is:"<<std::endl;
4282          int maxAcross=5;
4283          int limits[]={1,51,101,151,201,251,301,351,401};
4284          std::vector<std::string> types;
4285          types.push_back("Double parameters:");
4286          types.push_back("Branch and Cut double parameters:");
4287          types.push_back("Integer parameters:");
4288          types.push_back("Branch and Cut integer parameters:");
4289          types.push_back("Keyword parameters:");
4290          types.push_back("Branch and Cut keyword parameters:");
4291          types.push_back("Actions or string parameters:");
4292          types.push_back("Branch and Cut actions:");
4293          int iType;
4294          for (iType=0;iType<8;iType++) {
4295            int across=0;
4296            std::cout<<types[iType]<<"  ";
4297            for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4298              int type = parameters_[iParam].type();
4299              if (type>=limits[iType]
4300                  &&type<limits[iType+1]) {
4301                if (!across)
4302                  std::cout<<"  ";
4303                std::cout<<parameters_[iParam].matchName()<<"  ";
4304                across++;
4305                if (across==maxAcross) {
4306                  std::cout<<std::endl;
4307                  across=0;
4308                }
4309              }
4310            }
4311            if (across)
4312              std::cout<<std::endl;
4313          }
4314        } else if (type<101) {
4315          // get next field as double
4316          double value = CoinReadGetDoubleField(argc,argv,&valid);
4317          if (!valid) {
4318            if (type<51) {
4319              int returnCode;
4320              const char * message = 
4321                parameters_[iParam].setDoubleParameterWithMessage(lpSolver,value,returnCode);
4322              if (!noPrinting_&&strlen(message)) {
4323                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4324                  << message
4325                  <<CoinMessageEol;
4326              }
4327            } else if (type<81) {
4328              int returnCode;
4329              const char * message = 
4330                parameters_[iParam].setDoubleParameterWithMessage(model_,value,returnCode);
4331              if (!noPrinting_&&strlen(message)) {
4332                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4333                  << message
4334                  <<CoinMessageEol;
4335              }
4336            } else {
4337              int returnCode;
4338              const char * message = 
4339                parameters_[iParam].setDoubleParameterWithMessage(lpSolver,value,returnCode);
4340              if (!noPrinting_&&strlen(message)) {
4341                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4342                  << message
4343                  <<CoinMessageEol;
4344              }
4345              switch(type) {
4346              case DJFIX:
4347                djFix=value;
4348                if (goodModel&&djFix<1.0e20) {
4349                  // do some fixing
4350                  clpSolver = dynamic_cast< OsiClpSolverInterface*> (model_.solver());
4351                  clpSolver->initialSolve();
4352                  lpSolver = clpSolver->getModelPtr();
4353                  int numberColumns = lpSolver->numberColumns();
4354                  int i;
4355                  const char * type = lpSolver->integerInformation();
4356                  double * lower = lpSolver->columnLower();
4357                  double * upper = lpSolver->columnUpper();
4358                  double * solution = lpSolver->primalColumnSolution();
4359                  double * dj = lpSolver->dualColumnSolution();
4360                  int numberFixed=0;
4361                  double dextra4 = parameters_[whichParam(DEXTRA4,numberParameters_,parameters_)].doubleValue();
4362                  if (dextra4)
4363                    printf("Multiple for continuous dj fixing is %g\n",dextra4);
4364                  for (i=0;i<numberColumns;i++) {
4365                    double djValue = dj[i];
4366                    if (!type[i])
4367                      djValue *= dextra4;
4368                    if (type[i]||dextra4) {
4369                      double value = solution[i];
4370                      if (value<lower[i]+1.0e-5&&djValue>djFix) {
4371                        solution[i]=lower[i];
4372                        upper[i]=lower[i];
4373                        numberFixed++;
4374                      } else if (value>upper[i]-1.0e-5&&djValue<-djFix) {
4375                        solution[i]=upper[i];
4376                        lower[i]=upper[i];
4377                        numberFixed++;
4378                      }
4379                    }
4380                  }
4381                  sprintf(generalPrint,"%d columns fixed\n",numberFixed);
4382                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
4383                    << generalPrint
4384                    <<CoinMessageEol;
4385                }
4386                break;
4387              case GAPRATIO:
4388                gapRatio=value;
4389                break;
4390              case TIGHTENFACTOR:
4391                tightenFactor=value;
4392                if(!complicatedInteger)
4393                  defaultSettings=false; // user knows what she is doing
4394                break;
4395              default:
4396                break;
4397              }
4398            }
4399          } else if (valid==1) {
4400            std::cout<<" is illegal for double parameter "<<parameters_[iParam].name()<<" value remains "<<
4401              parameters_[iParam].doubleValue()<<std::endl;
4402          } else {
4403            std::cout<<parameters_[iParam].name()<<" has value "<<
4404              parameters_[iParam].doubleValue()<<std::endl;
4405          }
4406        } else if (type<201) {
4407          // get next field as int
4408          int value = CoinReadGetIntField(argc,argv,&valid);
4409          if (!valid) {
4410            if (type<151) {
4411              if (parameters_[iParam].type()==PRESOLVEPASS)
4412                preSolve = value;
4413              else if (parameters_[iParam].type()==IDIOT)
4414                doIdiot = value;
4415              else if (parameters_[iParam].type()==SPRINT)
4416                doSprint = value;
4417              else if (parameters_[iParam].type()==OUTPUTFORMAT)
4418                outputFormat = value;
4419              else if (parameters_[iParam].type()==SLPVALUE)
4420                slpValue = value;
4421              else if (parameters_[iParam].type()==CPP)
4422                cppValue = value;
4423              else if (parameters_[iParam].type()==PRESOLVEOPTIONS)
4424                presolveOptions = value;
4425              else if (parameters_[iParam].type()==PRINTOPTIONS)
4426                printOptions = value;
4427              else if (parameters_[iParam].type()==SUBSTITUTION)
4428                substitution = value;
4429              else if (parameters_[iParam].type()==DUALIZE)
4430                dualize = value;
4431              else if (parameters_[iParam].type()==PROCESSTUNE)
4432                tunePreProcess = value;
4433              else if (parameters_[iParam].type()==USESOLUTION)
4434                useSolution = value;
4435              else if (parameters_[iParam].type()==VERBOSE)
4436                verbose = value;
4437              int returnCode;
4438              const char * message = 
4439                parameters_[iParam].setIntParameterWithMessage(lpSolver,value,returnCode);
4440              if (!noPrinting_&&strlen(message)) {
4441                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4442                  << message
4443                  <<CoinMessageEol;
4444              }
4445            } else {
4446              if (parameters_[iParam].type()==CUTPASS)
4447                cutPass = value;
4448              else if (parameters_[iParam].type()==CUTPASSINTREE)
4449                cutPassInTree = value;
4450              else if (parameters_[iParam].type()==STRONGBRANCHING||
4451                       parameters_[iParam].type()==NUMBERBEFORE)
4452                strongChanged=true;
4453              else if (parameters_[iParam].type()==EXPERIMENT) {
4454                if (value>=1) {
4455                  if (!noPrinting_) {
4456                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4457                      <<"switching on dense factorization if small, and maybe fast fathoming"
4458                      <<CoinMessageEol;
4459                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4460                      <<"Gomory cuts using tolerance of 0.01 at root"
4461                      <<CoinMessageEol;
4462                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4463                      <<"Possible restart after 100 nodes if can fix many"
4464                      <<CoinMessageEol;
4465                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4466                      <<"extra options - -diving C -diveopt 3 -rins on -tune 6 -probing on -passf 30!"
4467                      <<CoinMessageEol;
4468                  }
4469                  // try changing tolerance at root
4470                  gomoryGen.setAwayAtRoot(0.01);
4471                  int iParam;
4472                  iParam = whichParam(DIVEOPT,numberParameters_,parameters_);
4473                  parameters_[iParam].setIntValue(3);
4474                  iParam = whichParam(FPUMPITS,numberParameters_,parameters_);
4475                  parameters_[iParam].setIntValue(30);
4476                  iParam = whichParam(PROCESSTUNE,numberParameters_,parameters_);
4477                  parameters_[iParam].setIntValue(6);
4478                  tunePreProcess=6;
4479                  iParam = whichParam(DIVINGC,numberParameters_,parameters_);
4480                  parameters_[iParam].setCurrentOption("on");
4481                  iParam = whichParam(RINS,numberParameters_,parameters_);
4482                  parameters_[iParam].setCurrentOption("on");
4483                  iParam = whichParam(PROBINGCUTS,numberParameters_,parameters_);
4484                  parameters_[iParam].setCurrentOption("on");
4485                  probingAction=1;
4486                }
4487              }
4488              int returnCode;
4489              const char * message = 
4490                parameters_[iParam].setIntParameterWithMessage(model_,value,returnCode);
4491              if (!noPrinting_&&strlen(message)) {
4492                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4493                  << message
4494                  <<CoinMessageEol;
4495              }
4496            }
4497          } else if (valid==1) {
4498            std::cout<<" is illegal for integer parameter "<<parameters_[iParam].name()<<" value remains "<<
4499              parameters_[iParam].intValue()<<std::endl;
4500          } else {
4501            std::cout<<parameters_[iParam].name()<<" has value "<<
4502              parameters_[iParam].intValue()<<std::endl;
4503          }
4504        } else if (type<301) {
4505          // one of several strings
4506          std::string value = CoinReadGetString(argc,argv);
4507          int action = parameters_[iParam].parameterOption(value);
4508          if (action<0) {
4509            if (value!="EOL") {
4510              // no match
4511              parameters_[iParam].printOptions();
4512            } else {
4513              // print current value
4514              std::cout<<parameters_[iParam].name()<<" has value "<<
4515                parameters_[iParam].currentOption()<<std::endl;
4516            }
4517          } else {
4518            const char * message = 
4519              parameters_[iParam].setCurrentOptionWithMessage(action);
4520            if (!noPrinting_&&strlen(message)) {
4521              generalMessageHandler->message(CLP_GENERAL,generalMessages)
4522                << message
4523                <<CoinMessageEol;
4524            }
4525            // for now hard wired
4526            switch (type) {
4527            case DIRECTION:
4528              if (action==0)
4529                lpSolver->setOptimizationDirection(1);
4530              else if (action==1)
4531                lpSolver->setOptimizationDirection(-1);
4532              else
4533                lpSolver->setOptimizationDirection(0);
4534              break;
4535            case DUALPIVOT:
4536              if (action==0) {
4537                ClpDualRowSteepest steep(3);
4538                lpSolver->setDualRowPivotAlgorithm(steep);
4539              } else if (action==1) {
4540                ClpDualRowDantzig dantzig;
4541                //ClpDualRowSteepest dantzig(5);
4542                lpSolver->setDualRowPivotAlgorithm(dantzig);
4543              } else if (action==2) {
4544                // partial steep
4545                ClpDualRowSteepest steep(2);
4546                lpSolver->setDualRowPivotAlgorithm(steep);
4547              } else {
4548                ClpDualRowSteepest steep;
4549                lpSolver->setDualRowPivotAlgorithm(steep);
4550              }
4551              break;
4552            case PRIMALPIVOT:
4553              if (action==0) {
4554                ClpPrimalColumnSteepest steep(3);
4555                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4556              } else if (action==1) {
4557                ClpPrimalColumnSteepest steep(0);
4558                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4559              } else if (action==2) {
4560                ClpPrimalColumnDantzig dantzig;
4561                lpSolver->setPrimalColumnPivotAlgorithm(dantzig);
4562              } else if (action==3) {
4563                ClpPrimalColumnSteepest steep(4);
4564                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4565              } else if (action==4) {
4566                ClpPrimalColumnSteepest steep(1);
4567                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4568              } else if (action==5) {
4569                ClpPrimalColumnSteepest steep(2);
4570                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4571              } else if (action==6) {
4572                ClpPrimalColumnSteepest steep(10);
4573                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4574              }
4575              break;
4576            case SCALING:
4577              lpSolver->scaling(action);
4578              solver->setHintParam(OsiDoScale,action!=0,OsiHintTry);
4579              doScaling = action;
4580              break;
4581            case AUTOSCALE:
4582              lpSolver->setAutomaticScaling(action!=0);
4583              break;
4584            case SPARSEFACTOR:
4585              lpSolver->setSparseFactorization((1-action)!=0);
4586              break;
4587            case BIASLU:
4588              lpSolver->factorization()->setBiasLU(action);
4589              break;
4590            case PERTURBATION:
4591              if (action==0)
4592                lpSolver->setPerturbation(50);
4593              else
4594                lpSolver->setPerturbation(100);
4595              break;
4596            case ERRORSALLOWED:
4597              allowImportErrors = action;
4598              break;
4599            case INTPRINT:
4600              printMode=action;
4601              break;
4602              //case ALGORITHM:
4603              //algorithm  = action;
4604              //defaultSettings=false; // user knows what she is doing
4605              //abort();
4606              //break;
4607            case KEEPNAMES:
4608              keepImportNames = 1-action;
4609              break;
4610            case PRESOLVE:
4611              if (action==0)
4612                preSolve = 5;
4613              else if (action==1)
4614                preSolve=0;
4615              else if (action==2)
4616                preSolve=10;
4617              else
4618                preSolveFile=true;
4619              break;
4620            case PFI:
4621              lpSolver->factorization()->setForrestTomlin(action==0);
4622              break;
4623            case CRASH:
4624              doCrash=action;
4625              break;
4626            case VECTOR:
4627              doVector=action;
4628              break;
4629            case MESSAGES:
4630              lpSolver->messageHandler()->setPrefix(action!=0);
4631              break;
4632            case CHOLESKY:
4633              choleskyType = action;
4634              break;
4635            case GAMMA:
4636              gamma=action;
4637              break;
4638            case BARRIERSCALE:
4639              scaleBarrier=action;
4640              break;
4641            case KKT:
4642              doKKT=action;
4643              break;
4644            case CROSSOVER:
4645              crossover=action;
4646              break;
4647            case SOS:
4648              doSOS=action;
4649              break;
4650            case GOMORYCUTS:
4651              defaultSettings=false; // user knows what she is doing
4652              gomoryAction = action;
4653              break;
4654            case PROBINGCUTS:
4655              defaultSettings=false; // user knows what she is doing
4656              probingAction = action;
4657              break;
4658            case KNAPSACKCUTS:
4659              defaultSettings=false; // user knows what she is doing
4660              knapsackAction = action;
4661              break;
4662            case REDSPLITCUTS:
4663              defaultSettings=false; // user knows what she is doing
4664              redsplitAction = action;
4665              break;
4666            case CLIQUECUTS:
4667              defaultSettings=false; // user knows what she is doing
4668              cliqueAction = action;
4669              break;
4670            case FLOWCUTS:
4671              defaultSettings=false; // user knows what she is doing
4672              flowAction = action;
4673              break;
4674            case MIXEDCUTS:
4675              defaultSettings=false; // user knows what she is doing
4676              mixedAction = action;
4677              break;
4678            case TWOMIRCUTS:
4679              defaultSettings=false; // user knows what she is doing
4680              twomirAction = action;
4681              break;
4682            case LANDPCUTS:
4683              defaultSettings=false; // user knows what she is doing
4684              landpAction = action;
4685              break;
4686            case RESIDCUTS:
4687              defaultSettings=false; // user knows what she is doing
4688              residualCapacityAction = action;
4689              break;
4690            case ROUNDING:
4691              defaultSettings=false; // user knows what she is doing
4692              break;
4693            case FPUMP:
4694              defaultSettings=false; // user knows what she is doing
4695              break;
4696            case RINS:
4697              break;
4698            case DINS:
4699              break;
4700            case RENS:
4701              break;
4702            case CUTSSTRATEGY:
4703              gomoryAction = action;
4704              probingAction = action;
4705              knapsackAction = action;
4706              cliqueAction = action;
4707              flowAction = action;
4708              mixedAction = action;
4709              twomirAction = action;
4710              //landpAction = action;
4711              parameters_[whichParam(GOMORYCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4712              parameters_[whichParam(PROBINGCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4713              parameters_[whichParam(KNAPSACKCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4714              parameters_[whichParam(CLIQUECUTS,numberParameters_,parameters_)].setCurrentOption(action);
4715              parameters_[whichParam(FLOWCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4716              parameters_[whichParam(MIXEDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4717              parameters_[whichParam(TWOMIRCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4718              if (!action) {
4719                redsplitAction = action;
4720                parameters_[whichParam(REDSPLITCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4721                landpAction = action;
4722                parameters_[whichParam(LANDPCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4723                residualCapacityAction = action;
4724                parameters_[whichParam(RESIDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4725              }
4726              break;
4727            case HEURISTICSTRATEGY:
4728              parameters_[whichParam(ROUNDING,numberParameters_,parameters_)].setCurrentOption(action);
4729              parameters_[whichParam(GREEDY,numberParameters_,parameters_)].setCurrentOption(action);
4730              parameters_[whichParam(COMBINE,numberParameters_,parameters_)].setCurrentOption(action);
4731              //parameters_[whichParam(LOCALTREE,numberParameters_,parameters_)].setCurrentOption(action);
4732              parameters_[whichParam(FPUMP,numberParameters_,parameters_)].setCurrentOption(action);
4733              break;
4734            case GREEDY:
4735            case DIVINGS:
4736            case DIVINGC:
4737            case DIVINGF:
4738            case DIVINGG:
4739            case DIVINGL:
4740            case DIVINGP:
4741            case DIVINGV:
4742            case COMBINE:
4743            case LOCALTREE:
4744              defaultSettings=false; // user knows what she is doing
4745              break;
4746            case COSTSTRATEGY:
4747              useCosts=action;
4748              break;
4749            case NODESTRATEGY:
4750              nodeStrategy=action;
4751              break;
4752            case PREPROCESS:
4753              preProcess = action;
4754              break;
4755              break;
4756            default:
4757              abort();
4758            }
4759          }
4760        } else {
4761          // action
4762          if (type==EXIT) {
4763#ifdef COIN_HAS_ASL
4764            if(statusUserFunction_[0]) {
4765              if (info.numberIntegers||info.numberBinary) {
4766                // integer
4767              } else {
4768                // linear
4769              }
4770              writeAmpl(&info);
4771              freeArrays2(&info);
4772              freeArgs(&info);
4773            }
4774#endif
4775            break; // stop all
4776          }
4777          switch (type) {
4778          case DUALSIMPLEX:
4779          case PRIMALSIMPLEX:
4780          case SOLVECONTINUOUS:
4781          case BARRIER:
4782            if (goodModel) {
4783              double objScale = 
4784                parameters_[whichParam(OBJSCALE2,numberParameters_,parameters_)].doubleValue();
4785              if (objScale!=1.0) {
4786                int iColumn;
4787                int numberColumns=lpSolver->numberColumns();
4788                double * dualColumnSolution = 
4789                  lpSolver->dualColumnSolution();
4790                ClpObjective * obj = lpSolver->objectiveAsObject();
4791                assert(dynamic_cast<ClpLinearObjective *> (obj));
4792                double offset;
4793                double * objective = obj->gradient(NULL,NULL,offset,true);
4794                for (iColumn=0;iColumn<numberColumns;iColumn++) {
4795                  dualColumnSolution[iColumn] *= objScale;
4796                  objective[iColumn] *= objScale;;
4797                }
4798                int iRow;
4799                int numberRows=lpSolver->numberRows();
4800                double * dualRowSolution = 
4801                  lpSolver->dualRowSolution();
4802                for (iRow=0;iRow<numberRows;iRow++) 
4803                  dualRowSolution[iRow] *= objScale;
4804                lpSolver->setObjectiveOffset(objScale*lpSolver->objectiveOffset());
4805              }
4806              ClpSolve::SolveType method;
4807              ClpSolve::PresolveType presolveType;
4808              ClpSimplex * model2 = lpSolver;
4809              if (dualize) {
4810                bool tryIt=true;
4811                double fractionColumn=1.0;
4812                double fractionRow=1.0;
4813                if (dualize==3) {
4814                  dualize=1;
4815                  int numberColumns=lpSolver->numberColumns();
4816                  int numberRows=lpSolver->numberRows();
4817                  if (numberRows<50000||5*numberColumns>numberRows) {
4818                    tryIt=false;
4819                  } else {
4820                    fractionColumn=0.1;
4821                    fractionRow=0.1;
4822                  }
4823                }
4824                if (tryIt) {
4825                  model2 = ((ClpSimplexOther *) model2)->dualOfModel(fractionRow,fractionColumn);
4826                  if (model2) {
4827                    sprintf(generalPrint,"Dual of model has %d rows and %d columns",
4828                            model2->numberRows(),model2->numberColumns());
4829                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4830                      << generalPrint
4831                      <<CoinMessageEol;
4832                    model2->setOptimizationDirection(1.0);
4833                  } else {
4834                    model2 = lpSolver;
4835                    dualize=0;
4836                  }
4837                } else {
4838                  dualize=0;
4839                }
4840              }
4841              if (noPrinting_)
4842                lpSolver->setLogLevel(0);
4843              ClpSolve solveOptions;
4844              solveOptions.setPresolveActions(presolveOptions);
4845              solveOptions.setSubstitution(substitution);
4846              if (preSolve!=5&&preSolve) {
4847                presolveType=ClpSolve::presolveNumber;
4848                if (preSolve<0) {
4849                  preSolve = - preSolve;
4850                  if (preSolve<=100) {
4851                    presolveType=ClpSolve::presolveNumber;
4852                    sprintf(generalPrint,"Doing %d presolve passes - picking up non-costed slacks",
4853                           preSolve);
4854                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4855                      << generalPrint
4856                      <<CoinMessageEol;
4857                    solveOptions.setDoSingletonColumn(true);
4858                  } else {
4859                    preSolve -=100;
4860                    presolveType=ClpSolve::presolveNumberCost;
4861                    sprintf(generalPrint,"Doing %d presolve passes - picking up costed slacks",
4862                           preSolve);
4863                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4864                      << generalPrint
4865                      <<CoinMessageEol;
4866                  }
4867                } 
4868              } else if (preSolve) {
4869                presolveType=ClpSolve::presolveOn;
4870              } else {
4871                presolveType=ClpSolve::presolveOff;
4872              }
4873              solveOptions.setPresolveType(presolveType,preSolve);
4874              if (type==DUALSIMPLEX||type==SOLVECONTINUOUS) {
4875                method=ClpSolve::useDual;
4876              } else if (type==PRIMALSIMPLEX) {
4877                method=ClpSolve::usePrimalorSprint;
4878              } else {
4879                method = ClpSolve::useBarrier;
4880                if (crossover==1) {
4881                  method=ClpSolve::useBarrierNoCross;
4882                } else if (crossover==2) {
4883                  ClpObjective * obj = lpSolver->objectiveAsObject();
4884                  if (obj->type()>1) {
4885                    method=ClpSolve::useBarrierNoCross;
4886                    presolveType=ClpSolve::presolveOff;
4887                    solveOptions.setPresolveType(presolveType,preSolve);
4888                  } 
4889                }
4890              }
4891              solveOptions.setSolveType(method);
4892              if(preSolveFile)
4893                presolveOptions |= 0x40000000;
4894              solveOptions.setSpecialOption(4,presolveOptions);
4895              solveOptions.setSpecialOption(5,printOptions);
4896              if (doVector) {
4897                ClpMatrixBase * matrix = lpSolver->clpMatrix();
4898                if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
4899                  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
4900                  clpMatrix->makeSpecialColumnCopy();
4901                }
4902              }
4903              if (method==ClpSolve::useDual) {
4904                // dual
4905                if (doCrash)
4906                  solveOptions.setSpecialOption(0,1,doCrash); // crash
4907                else if (doIdiot)
4908                  solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
4909              } else if (method==ClpSolve::usePrimalorSprint) {
4910                // primal
4911                // if slp turn everything off
4912                if (slpValue>0) {
4913                  doCrash=false;
4914                  doSprint=0;
4915                  doIdiot=-1;
4916                  solveOptions.setSpecialOption(1,10,slpValue); // slp
4917                  method=ClpSolve::usePrimal;
4918                }
4919                if (doCrash) {
4920                  solveOptions.setSpecialOption(1,1,doCrash); // crash
4921                } else if (doSprint>0) {
4922                  // sprint overrides idiot
4923                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
4924                } else if (doIdiot>0) {
4925                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
4926                } else if (slpValue<=0) {
4927                  if (doIdiot==0) {
4928                    if (doSprint==0)
4929                      solveOptions.setSpecialOption(1,4); // all slack
4930                    else
4931                      solveOptions.setSpecialOption(1,9); // all slack or sprint
4932                  } else {
4933                    if (doSprint==0)
4934                      solveOptions.setSpecialOption(1,8); // all slack or idiot
4935                    else
4936                      solveOptions.setSpecialOption(1,7); // initiative
4937                  }
4938                }
4939                if (basisHasValues==-1)
4940                  solveOptions.setSpecialOption(1,11); // switch off values
4941              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
4942                int barrierOptions = choleskyType;
4943                if (scaleBarrier)
4944                  barrierOptions |= 8;
4945                if (doKKT)
4946                  barrierOptions |= 16;
4947                if (gamma)
4948                  barrierOptions |= 32*gamma;
4949                if (crossover==3) 
4950                  barrierOptions |= 256; // try presolve in crossover
4951                solveOptions.setSpecialOption(4,barrierOptions);
4952              }
4953              model2->setMaximumSeconds(model_.getMaximumSeconds());
4954#ifdef COIN_HAS_LINK
4955              OsiSolverInterface * coinSolver = model_.solver();
4956              OsiSolverLink * linkSolver = dynamic_cast< OsiSolverLink*> (coinSolver);
4957              if (!linkSolver) {
4958                model2->initialSolve(solveOptions);
4959              } else {
4960                // special solver
4961                int testOsiOptions = parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].intValue();
4962                double * solution = NULL;
4963                if (testOsiOptions<10) {
4964                  solution = linkSolver->nonlinearSLP(slpValue>0 ? slpValue : 20 ,1.0e-5);
4965                } else if (testOsiOptions>=10) {
4966                  CoinModel coinModel = *linkSolver->coinModel();
4967                  ClpSimplex * tempModel = approximateSolution(coinModel,slpValue>0 ? slpValue : 50 ,1.0e-5,0);
4968                  assert (tempModel);
4969                  solution = CoinCopyOfArray(tempModel->primalColumnSolution(),coinModel.numberColumns());
4970                  model2->setObjectiveValue(tempModel->objectiveValue());
4971                  model2->setProblemStatus(tempModel->problemStatus());
4972                  model2->setSecondaryStatus(tempModel->secondaryStatus());
4973                  delete tempModel;
4974                }
4975                if (solution) {
4976                  memcpy(model2->primalColumnSolution(),solution,
4977                         CoinMin(model2->numberColumns(),linkSolver->coinModel()->numberColumns())*sizeof(double));
4978                  delete [] solution;
4979                } else {
4980                  printf("No nonlinear solution\n");
4981                }
4982              }
4983#else
4984              model2->initialSolve(solveOptions);
4985#endif
4986              {
4987                // map states
4988                /* clp status
4989                   -1 - unknown e.g. before solve or if postSolve says not optimal
4990                   0 - optimal
4991                   1 - primal infeasible
4992                   2 - dual infeasible
4993                   3 - stopped on iterations or time
4994                   4 - stopped due to errors
4995                   5 - stopped by event handler (virtual int ClpEventHandler::event()) */
4996                /* cbc status
4997                   -1 before branchAndBound
4998                   0 finished - check isProvenOptimal or isProvenInfeasible to see if solution found
4999                   (or check value of best solution)
5000                   1 stopped - on maxnodes, maxsols, maxtime
5001                   2 difficulties so run was abandoned
5002                   (5 event user programmed event occurred) */
5003                /* clp secondary status of problem - may get extended
5004                   0 - none
5005                   1 - primal infeasible because dual limit reached OR probably primal
5006                   infeasible but can't prove it (main status 4)
5007                   2 - scaled problem optimal - unscaled problem has primal infeasibilities
5008                   3 - scaled problem optimal - unscaled problem has dual infeasibilities
5009                   4 - scaled problem optimal - unscaled problem has primal and dual infeasibilities
5010                   5 - giving up in primal with flagged variables
5011                   6 - failed due to empty problem check
5012                   7 - postSolve says not optimal
5013                   8 - failed due to bad element check
5014                   9 - status was 3 and stopped on time
5015                   100 up - translation of enum from ClpEventHandler
5016                */
5017                /* cbc secondary status of problem
5018                   -1 unset (status_ will also be -1)
5019                   0 search completed with solution
5020                   1 linear relaxation not feasible (or worse than cutoff)
5021                   2 stopped on gap
5022                   3 stopped on nodes
5023                   4 stopped on time
5024                   5 stopped on user event
5025                   6 stopped on solutions
5026                   7 linear relaxation unbounded
5027                */
5028                int iStatus = model2->status();
5029                int iStatus2 = model2->secondaryStatus();
5030                if (iStatus==0) {
5031                  iStatus2=0;
5032                } else if (iStatus==1) {
5033                  iStatus=0;
5034                  iStatus2=1; // say infeasible
5035                } else if (iStatus==2) {
5036                  iStatus=0;
5037                  iStatus2=7; // say unbounded
5038                } else if (iStatus==3) {
5039                  iStatus=1;
5040                  if (iStatus2==9)
5041                    iStatus2=4;
5042                  else
5043                    iStatus2=3; // Use nodes - as closer than solutions
5044                } else if (iStatus==4) {
5045                  iStatus=2; // difficulties
5046                  iStatus2=0; 
5047                }
5048                model_.setProblemStatus(iStatus);
5049                model_.setSecondaryStatus(iStatus2);
5050                assert (lpSolver==clpSolver->getModelPtr());
5051                assert (clpSolver==model_.solver());
5052                clpSolver->setWarmStart(NULL);
5053                // and in babModel if exists
5054                if (babModel_) {
5055                  babModel_->setProblemStatus(iStatus);
5056                  babModel_->setSecondaryStatus(iStatus2);
5057                } 
5058#if NEW_STYLE_SOLVER
5059                int returnCode = callBack_->callBack(&model_,1);
5060#else
5061                int returnCode=callBack(&model,1);
5062#endif
5063                if (returnCode) {
5064                  // exit if user wants
5065#if NEW_STYLE_SOLVER
5066                  updateModel(model2,returnMode);
5067#else
5068                  delete babModel_;
5069                  babModel_ = NULL;
5070#endif
5071                  return returnCode;
5072                }
5073              }
5074              basisHasValues=1;
5075              if (dualize) {
5076                int returnCode=((ClpSimplexOther *) lpSolver)->restoreFromDual(model2);
5077                if (model2->status()==3)
5078                  returnCode=0;
5079                delete model2;
5080                if (returnCode&&dualize!=2)
5081                  lpSolver->primal(1);
5082                model2=lpSolver;
5083              }
5084#if NEW_STYLE_SOLVER
5085              updateModel(model2,returnMode);
5086              for (iUser=0;iUser<numberUserFunctions_;iUser++) {
5087                if (statusUserFunction_[iUser])
5088                  userFunction_[iUser]->exportSolution(this,1);
5089              }
5090#endif
5091#ifdef COIN_HAS_ASL
5092              if (statusUserFunction_[0]) {
5093                double value = model2->getObjValue()*model2->getObjSense();
5094                char buf[300];
5095                int pos=0;
5096                int iStat = model2->status();
5097                if (iStat==0) {
5098                  pos += sprintf(buf+pos,"optimal," );
5099                } else if (iStat==1) {
5100                  // infeasible
5101                  pos += sprintf(buf+pos,"infeasible,");
5102                } else if (iStat==2) {
5103                  // unbounded
5104                  pos += sprintf(buf+pos,"unbounded,");
5105                } else if (iStat==3) {
5106                  pos += sprintf(buf+pos,"stopped on iterations or time,");
5107                } else if (iStat==4) {
5108                  iStat = 7;
5109                  pos += sprintf(buf+pos,"stopped on difficulties,");
5110                } else if (iStat==5) {
5111                  iStat = 3;
5112                  pos += sprintf(buf+pos,"stopped on ctrl-c,");
5113                } else {
5114                  pos += sprintf(buf+pos,"status unknown,");
5115                  iStat=6;
5116                }
5117                info.problemStatus=iStat;
5118                info.objValue = value;
5119                pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
5120                               value);
5121                sprintf(buf+pos,"\n%d iterations",
5122                        model2->getIterationCount());
5123                free(info.primalSolution);
5124                int numberColumns=model2->numberColumns();
5125                info.primalSolution = (double *) malloc(numberColumns*sizeof(double));
5126                CoinCopyN(model2->primalColumnSolution(),numberColumns,info.primalSolution);
5127                int numberRows = model2->numberRows();
5128                free(info.dualSolution);
5129                info.dualSolution = (double *) malloc(numberRows*sizeof(double));
5130                CoinCopyN(model2->dualRowSolution(),numberRows,info.dualSolution);
5131                CoinWarmStartBasis * basis = model2->getBasis();
5132                free(info.rowStatus);
5133                info.rowStatus = (int *) malloc(numberRows*sizeof(int));
5134                free(info.columnStatus);
5135                info.columnStatus = (int *) malloc(numberColumns*sizeof(int));
5136                // Put basis in
5137                int i;
5138                // free,basic,ub,lb are 0,1,2,3
5139                for (i=0;i<numberRows;i++) {
5140                  CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
5141                  info.rowStatus[i]=status;
5142                }
5143                for (i=0;i<numberColumns;i++) {
5144                  CoinWarmStartBasis::Status status = basis->getStructStatus(i);
5145                  info.columnStatus[i]=status;
5146                }
5147                // put buffer into info
5148                strcpy(info.buffer,buf);
5149                delete basis;
5150              }
5151#endif
5152            } else {
5153#ifndef DISALLOW_PRINTING
5154              std::cout<<"** Current model not valid"<<std::endl;
5155#endif
5156            }
5157            break;
5158          case STATISTICS:
5159            if (goodModel) {
5160              // If presolve on look at presolved
5161              bool deleteModel2=false;
5162              ClpSimplex * model2 = lpSolver;
5163              if (preSolve) {
5164                ClpPresolve pinfo;
5165                int presolveOptions2 = presolveOptions&~0x40000000;
5166                if ((presolveOptions2&0xffff)!=0)
5167                  pinfo.setPresolveActions(presolveOptions2);
5168                pinfo.setSubstitution(substitution);
5169                if ((printOptions&1)!=0)
5170                  pinfo.statistics();
5171                double presolveTolerance = 
5172                  parameters_[whichParam(PRESOLVETOLERANCE,numberParameters_,parameters_)].doubleValue();
5173                model2 = 
5174                  pinfo.presolvedModel(*lpSolver,presolveTolerance,
5175                                       true,preSolve);
5176                if (model2) {
5177                  printf("Statistics for presolved model\n");
5178                  deleteModel2=true;
5179                } else {
5180                  printf("Presolved model looks infeasible - will use unpresolved\n");
5181                  model2 = lpSolver;
5182                }
5183              } else {
5184                printf("Statistics for unpresolved model\n");
5185                model2 =  lpSolver;
5186              }
5187              statistics(lpSolver,model2);
5188              if (deleteModel2)
5189                delete model2;
5190            } else {
5191#ifndef DISALLOW_PRINTING
5192              std::cout<<"** Current model not valid"<<std::endl;
5193#endif
5194            }
5195            break;
5196          case TIGHTEN:
5197            if (goodModel) {
5198              int numberInfeasibilities = lpSolver->tightenPrimalBounds();
5199              if (numberInfeasibilities)
5200                std::cout<<"** Analysis indicates model infeasible"<<std::endl;
5201            } else {
5202#ifndef DISALLOW_PRINTING
5203              std::cout<<"** Current model not valid"<<std::endl;
5204#endif
5205            }
5206            break;
5207          case PLUSMINUS:
5208            if (goodModel) {
5209              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
5210              ClpPackedMatrix* clpMatrix =
5211                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
5212              if (clpMatrix) {
5213                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
5214                if (newMatrix->getIndices()) {
5215                  lpSolver->replaceMatrix(newMatrix);
5216                  delete saveMatrix;
5217                  std::cout<<"Matrix converted to +- one matrix"<<std::endl;
5218                } else {
5219                  std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
5220                }
5221              } else {
5222                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
5223              }
5224            } else {
5225#ifndef DISALLOW_PRINTING
5226              std::cout<<"** Current model not valid"<<std::endl;
5227#endif
5228            }
5229            break;
5230          case OUTDUPROWS:
5231            if (goodModel) {
5232              int numberRows = clpSolver->getNumRows();
5233              //int nOut = outDupRow(clpSolver);
5234              CglDuplicateRow dupcuts(clpSolver);
5235              storedCuts = dupcuts.outDuplicates(clpSolver)!=0;
5236              int nOut = numberRows-clpSolver->getNumRows();
5237              if (nOut&&!noPrinting_)
5238                sprintf(generalPrint,"%d rows eliminated",nOut);
5239              generalMessageHandler->message(CLP_GENERAL,generalMessages)
5240                << generalPrint
5241                <<CoinMessageEol;
5242            } else {
5243#ifndef DISALLOW_PRINTING
5244              std::cout<<"** Current model not valid"<<std::endl;
5245#endif
5246            }
5247            break;
5248          case NETWORK:
5249            if (goodModel) {
5250              ClpMatrixBase * saveMatrix = lpSolver->clpMatrix();
5251              ClpPackedMatrix* clpMatrix =
5252                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
5253              if (clpMatrix) {
5254                ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
5255                if (newMatrix->getIndices()) {
5256                  lpSolver->replaceMatrix(newMatrix);
5257                  delete saveMatrix;
5258                  std::cout<<"Matrix converted to network matrix"<<std::endl;
5259                } else {
5260                  std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
5261                }
5262              } else {
5263                std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
5264              }
5265            } else {
5266#ifndef DISALLOW_PRINTING
5267              std::cout<<"** Current model not valid"<<std::endl;
5268#endif
5269            }
5270            break;
5271          case DOHEURISTIC:
5272            if (goodModel) {
5273              int vubAction = parameters_[whichParam(VUBTRY,numberParameters_,parameters_)].intValue();
5274              if (vubAction!=-1) {
5275                // look at vubs
5276                // Just ones which affect >= extra3
5277                int extra3 = parameters_[whichParam(EXTRA3,numberParameters_,parameters_)].intValue();
5278                /* 2 is cost above which to fix if feasible
5279                   3 is fraction of integer variables fixed (0.97)
5280                   4 is fraction of all variables fixed (0.0)
5281                */
5282                double dextra[5];
5283                dextra[1] = parameters_[whichParam(DEXTRA1,numberParameters_,parameters_)].doubleValue();
5284                dextra[2] = parameters_[whichParam(DEXTRA2,numberParameters_,parameters_)].doubleValue();
5285                dextra[3] = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
5286                dextra[4] = parameters_[whichParam(DEXTRA4,numberParameters_,parameters_)].doubleValue();
5287                if (!dextra[3])
5288                  dextra[3] = 0.97;
5289                //OsiClpSolverInterface * newSolver =
5290                fixVubs(model_,extra3,vubAction,generalMessageHandler,
5291                        debugValues,dextra);
5292                //assert (!newSolver);
5293              }
5294              // Actually do heuristics
5295              doHeuristics(&model_,2);
5296              if (model_.bestSolution()) {
5297                model_.setProblemStatus(1);
5298                model_.setSecondaryStatus(6);
5299#ifdef COIN_HAS_ASL
5300                if (statusUserFunction_[0]) {
5301                  double value = model_.getObjValue();
5302                  char buf[300];
5303                  int pos=0;
5304                  pos += sprintf(buf+pos,"feasible,");
5305                  info