source: stable/2.0/Cbc/src/CbcSolver.cpp @ 806

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

make sure says 2.0

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