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

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

fixes for bonmin and other stuff

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