source: stable/2.1/Cbc/src/CbcSolver.cpp @ 957

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

fix few bugs

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