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

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

changes

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