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

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

changes to fpump

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