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

Last change on this file since 824 was 824, checked in by forrest, 12 years ago

for unit test return code

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