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

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

add ampl stuff after heuristic

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