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

Last change on this file since 1132 was 1132, checked in by forrest, 10 years ago

chnages to try and make faster

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