source: branches/heur/Cbc/src/CbcSolver.cpp @ 885

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

modify csv to add cuts

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