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

Last change on this file since 822 was 822, checked in by forrest, 13 years ago

prepare to use cliques

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