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

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

changes to heuristics and allow collection of more statistics

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