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

Last change on this file since 1024 was 1024, checked in by ladanyi, 11 years ago

fixed externals and removed references to the DINS heur as the source of that is not yet checked in

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