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

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

improve hot start

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