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

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

add diving heuristics

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