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

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

try and make faster

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