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

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

for ratiogap

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