source: stable/2.3/Cbc/src/CbcSolver.cpp @ 1234

Last change on this file since 1234 was 1234, checked in by forrest, 11 years ago

say infeasible in CbcSolver? output

File size: 377.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 not in integer
3773  int integerStatus=-1;
3774  // Say no resolve after cuts
3775  model_.setResolveAfterTakeOffCuts(false);
3776  // see if log in list
3777  for (int i=1;i<argc;i++) {
3778    if (!strncmp(argv[i],"log",3)) {
3779      const char * equals = strchr(argv[i],'=');
3780      if (equals&&atoi(equals+1)!=0) 
3781        noPrinting_=false;
3782      else
3783        noPrinting_=true;
3784      break;
3785    } else if (!strncmp(argv[i],"-log",4)&&i<argc-1) {
3786      if (atoi(argv[i+1])!=0) 
3787        noPrinting_=false;
3788      else
3789        noPrinting_=true;
3790      break;
3791    }
3792  }
3793  double time0;
3794  {
3795    double time1 = CoinCpuTime(),time2;
3796    time0=time1;
3797    bool goodModel=(originalSolver->getNumCols()) ? true : false;
3798
3799    // register signal handler
3800    //CoinSighandler_t saveSignal=signal(SIGINT,signal_handler);
3801    signal(SIGINT,signal_handler);
3802    // Set up all non-standard stuff
3803    int cutPass=-1234567;
3804    int cutPassInTree=-1234567;
3805    int tunePreProcess=0;
3806    int testOsiParameters=-1;
3807    // 0 normal, 1 from ampl or MIQP etc (2 allows cuts)
3808    int complicatedInteger=0;
3809    OsiSolverInterface * solver = model_.solver();
3810    if (noPrinting_) 
3811      setCbcOrClpPrinting(false);
3812#ifndef CBC_OTHER_SOLVER
3813    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3814    ClpSimplex * lpSolver = clpSolver->getModelPtr();
3815    if (noPrinting_) {
3816      lpSolver->setLogLevel(0);
3817    }
3818#else
3819    ClpSimplex * lpSolver = NULL;
3820#endif
3821    // For priorities etc
3822    int * priorities=NULL;
3823    int * branchDirection=NULL;
3824    double * pseudoDown=NULL;
3825    double * pseudoUp=NULL;
3826    double * solutionIn = NULL;
3827    int * prioritiesIn = NULL;
3828    int numberSOS = 0;
3829    int * sosStart = NULL;
3830    int * sosIndices = NULL;
3831    char * sosType = NULL;
3832    double * sosReference = NULL;
3833    int * cut=NULL;
3834    int * sosPriority=NULL;
3835    CglStored storedAmpl;
3836    CoinModel * coinModel = NULL;
3837    CoinModel saveCoinModel;
3838    CoinModel saveTightenedModel;
3839    int * whichColumn = NULL;
3840    int * knapsackStart=NULL;
3841    int * knapsackRow=NULL;
3842    int numberKnapsack=0;
3843#if NEW_STYLE_SOLVER
3844    int numberInputs=0;
3845    readMode_=CbcOrClpRead_mode;
3846    for (iUser=0;iUser<numberUserFunctions_;iUser++) {
3847      int status = userFunction_[iUser]->importData(this,argc,const_cast<char **>(argv));
3848      if (status>=0) {
3849        if (!status) {
3850          numberInputs++;
3851          statusUserFunction_[iUser]=1;
3852          goodModel=true;
3853          solver = model_.solver();
3854#ifndef CBC_OTHER_SOLVER
3855          clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3856          lpSolver = clpSolver->getModelPtr();
3857#endif
3858        } else {
3859          printf("Bad input from user function %s\n",userFunction_[iUser]->name().c_str());
3860          abort();
3861        }
3862      }
3863    }
3864    if (numberInputs>1) {
3865      printf("Two or more user inputs!\n");
3866      abort();
3867    }
3868    if (originalCoinModel_)
3869      complicatedInteger=1;
3870    testOsiParameters = intValue(TESTOSI);
3871    if (noPrinting_) {
3872      model_.messageHandler()->setLogLevel(0);
3873      setCbcOrClpPrinting(false);
3874    }
3875    CbcOrClpRead_mode=readMode_;
3876#endif
3877#ifdef COIN_HAS_ASL
3878    ampl_info info;
3879    {
3880    memset(&info,0,sizeof(info));
3881    if (argc>2&&!strcmp(argv[2],"-AMPL")) {
3882      statusUserFunction_[0]=1;
3883      // see if log in list
3884      noPrinting_=true;
3885      for (int i=1;i<argc;i++) {
3886        if (!strncmp(argv[i],"log",3)) {
3887          const char * equals = strchr(argv[i],'=');
3888          if (equals&&atoi(equals+1)>0) {
3889            noPrinting_=false;
3890            info.logLevel=atoi(equals+1);
3891            int log = whichParam(LOGLEVEL,numberParameters_,parameters_);
3892            parameters_[log].setIntValue(info.logLevel);
3893            // mark so won't be overWritten
3894            info.numberRows=-1234567;
3895            break;
3896          }
3897        }
3898      }
3899
3900      union { void * voidModel; CoinModel * model; } coinModelStart;
3901      coinModelStart.model=NULL;
3902      int returnCode = readAmpl(&info,argc, const_cast<char **>(argv),& coinModelStart.voidModel);
3903      coinModel=coinModelStart.model;
3904      if (returnCode)
3905        return returnCode;
3906      CbcOrClpRead_mode=2; // so will start with parameters
3907      // see if log in list (including environment)
3908      for (int i=1;i<info.numberArguments;i++) {
3909        if (!strcmp(info.arguments[i],"log")) {
3910          if (i<info.numberArguments-1&&atoi(info.arguments[i+1])>0)
3911            noPrinting_=false;
3912          break;
3913        }
3914      }
3915      if (noPrinting_) {
3916        model_.messageHandler()->setLogLevel(0);
3917        setCbcOrClpPrinting(false);
3918      }
3919      if (!noPrinting_)
3920        printf("%d rows, %d columns and %d elements\n",
3921               info.numberRows,info.numberColumns,info.numberElements);
3922#ifdef COIN_HAS_LINK
3923      if (!coinModel) {
3924#endif
3925      solver->loadProblem(info.numberColumns,info.numberRows,info.starts,
3926                          info.rows,info.elements,
3927                          info.columnLower,info.columnUpper,info.objective,
3928                          info.rowLower,info.rowUpper);
3929      // take off cuts if ampl wants that
3930      if (info.cut&&0) {
3931        printf("AMPL CUTS OFF until global cuts fixed\n");
3932        info.cut=NULL;
3933      }
3934      if (info.cut) {
3935        int numberRows = info.numberRows;
3936        int * whichRow = new int [numberRows];
3937        // Row copy
3938        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
3939        const double * elementByRow = matrixByRow->getElements();
3940        const int * column = matrixByRow->getIndices();
3941        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
3942        const int * rowLength = matrixByRow->getVectorLengths();
3943       
3944        const double * rowLower = solver->getRowLower();
3945        const double * rowUpper = solver->getRowUpper();
3946        int nDelete=0;
3947        for (int iRow=0;iRow<numberRows;iRow++) {
3948          if (info.cut[iRow]) {
3949            whichRow[nDelete++]=iRow;
3950            int start = rowStart[iRow];
3951            storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
3952                          rowLength[iRow],column+start,elementByRow+start);
3953          }
3954        }
3955        solver->deleteRows(nDelete,whichRow);
3956        delete [] whichRow;
3957      }
3958#ifdef COIN_HAS_LINK
3959      } else {
3960#ifndef CBC_OTHER_SOLVER
3961        // save
3962        saveCoinModel = *coinModel;
3963        // load from coin model
3964        OsiSolverLink solver1;
3965        OsiSolverInterface * solver2 = solver1.clone();
3966        model_.assignSolver(solver2,false);
3967        OsiSolverLink * si =
3968          dynamic_cast<OsiSolverLink *>(model_.solver()) ;
3969        assert (si != NULL);
3970        si->setDefaultMeshSize(0.001);
3971        // need some relative granularity
3972        si->setDefaultBound(100.0);
3973        double dextra3 = parameters_[whichParam(DEXTRA3,numberParameters_,parameters_)].doubleValue();
3974        if (dextra3)
3975          si->setDefaultMeshSize(dextra3);
3976        si->setDefaultBound(100000.0);
3977        si->setIntegerPriority(1000);
3978        si->setBiLinearPriority(10000);
3979        CoinModel * model2 = reinterpret_cast<CoinModel *> (coinModel);
3980        int logLevel = parameters_[whichParam(LOGLEVEL,numberParameters_,parameters_)].intValue();
3981        si->load(*model2,true,logLevel);
3982        // redo
3983        solver = model_.solver();
3984        clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
3985        lpSolver = clpSolver->getModelPtr();
3986        clpSolver->messageHandler()->setLogLevel(0) ;
3987        testOsiParameters=0;
3988        parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].setIntValue(0);
3989        complicatedInteger=1;
3990        if (info.cut) {
3991          printf("Sorry - can't do cuts with LOS as ruins delicate row order\n");
3992          abort();
3993          int numberRows = info.numberRows;
3994          int * whichRow = new int [numberRows];
3995          // Row copy
3996          const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
3997          const double * elementByRow = matrixByRow->getElements();
3998          const int * column = matrixByRow->getIndices();
3999          const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
4000          const int * rowLength = matrixByRow->getVectorLengths();
4001         
4002          const double * rowLower = solver->getRowLower();
4003          const double * rowUpper = solver->getRowUpper();
4004          int nDelete=0;
4005          for (int iRow=0;iRow<numberRows;iRow++) {
4006            if (info.cut[iRow]) {
4007              whichRow[nDelete++]=iRow;
4008              int start = rowStart[iRow];
4009              storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
4010                                rowLength[iRow],column+start,elementByRow+start);
4011            }
4012          }
4013          solver->deleteRows(nDelete,whichRow);
4014          // and special matrix
4015          si->cleanMatrix()->deleteRows(nDelete,whichRow);
4016          delete [] whichRow;
4017        }
4018#endif
4019      }
4020#endif
4021      // If we had a solution use it
4022      if (info.primalSolution) {
4023        solver->setColSolution(info.primalSolution);
4024      }
4025      // status
4026      if (info.rowStatus) {
4027        unsigned char * statusArray = lpSolver->statusArray();
4028        int i;
4029        for (i=0;i<info.numberColumns;i++)
4030          statusArray[i]= info.columnStatus[i];
4031        statusArray+=info.numberColumns;
4032        for (i=0;i<info.numberRows;i++)
4033          statusArray[i]= info.rowStatus[i];
4034        CoinWarmStartBasis * basis = lpSolver->getBasis();
4035        solver->setWarmStart(basis);
4036        delete basis;
4037      }
4038      freeArrays1(&info);
4039      // modify objective if necessary
4040      solver->setObjSense(info.direction);
4041      solver->setDblParam(OsiObjOffset,info.offset);
4042      if (info.offset) {
4043        sprintf(generalPrint,"Ampl objective offset is %g",
4044                info.offset);
4045        generalMessageHandler->message(CLP_GENERAL,generalMessages)
4046          << generalPrint
4047          <<CoinMessageEol;
4048      }
4049      // Set integer variables (unless nonlinear when set)
4050      if (!info.nonLinear) {
4051        for (int i=info.numberColumns-info.numberIntegers;
4052             i<info.numberColumns;i++)
4053          solver->setInteger(i);
4054      }
4055      goodModel=true;
4056      // change argc etc
4057      argc = info.numberArguments;
4058      argv = const_cast<const char **>(info.arguments);
4059    }
4060    }
4061#endif   
4062    // default action on import
4063    int allowImportErrors=0;
4064    int keepImportNames=1;
4065    int doIdiot=-1;
4066    int outputFormat=2;
4067    int slpValue=-1;
4068    int cppValue=-1;
4069    int printOptions=0;
4070    int printMode=0;
4071    int presolveOptions=0;
4072    int substitution=3;
4073    int dualize=3;
4074    int doCrash=0;
4075    int doVector=0;
4076    int doSprint=-1;
4077    int doScaling=4;
4078    // set reasonable defaults
4079    int preSolve=5;
4080    int preProcess=4;
4081    bool useStrategy=false;
4082    bool preSolveFile=false;
4083    bool strongChanged=false;
4084    bool pumpChanged=false;
4085   
4086    double djFix=1.0e100;
4087    double tightenFactor=0.0;
4088    const char dirsep =  CoinFindDirSeparator();
4089    std::string directory;
4090    std::string dirSample;
4091    std::string dirNetlib;
4092    std::string dirMiplib;
4093    if (dirsep == '/') {
4094      directory = "./";
4095      dirSample = "../../Data/Sample/";
4096      dirNetlib = "../../Data/Netlib/";
4097      dirMiplib = "../../Data/miplib3/";
4098    } else {
4099      directory = ".\\";
4100      dirSample = "..\\..\\Data\\Sample\\";
4101      dirNetlib = "..\\..\\Data\\Netlib\\";
4102      dirMiplib = "..\\..\\Data\\miplib3\\";
4103    }
4104    std::string defaultDirectory = directory;
4105    std::string importFile ="";
4106    std::string exportFile ="default.mps";
4107    std::string importBasisFile ="";
4108    std::string importPriorityFile ="";
4109    std::string debugFile="";
4110    std::string printMask="";
4111    double * debugValues = NULL;
4112    int numberDebugValues = -1;
4113    int basisHasValues=0;
4114    std::string exportBasisFile ="default.bas";
4115    std::string saveFile ="default.prob";
4116    std::string restoreFile ="default.prob";
4117    std::string solutionFile ="stdout";
4118    std::string solutionSaveFile ="solution.file";
4119    int slog = whichParam(SOLVERLOGLEVEL,numberParameters_,parameters_);
4120    int log = whichParam(LOGLEVEL,numberParameters_,parameters_);
4121#ifndef CBC_OTHER_SOLVER
4122    double normalIncrement=model_.getCutoffIncrement();;
4123#endif
4124    if (testOsiParameters>=0) {
4125      // trying nonlinear - switch off some stuff
4126      preProcess=0;
4127    }
4128    // Set up likely cut generators and defaults
4129    int nodeStrategy=0;
4130    int doSOS=1;
4131    int verbose=0;
4132    CglGomory gomoryGen;
4133    // try larger limit
4134    gomoryGen.setLimitAtRoot(1000);
4135    gomoryGen.setLimit(50);
4136    // set default action (0=off,1=on,2=root)
4137    int gomoryAction=3;
4138
4139    CglProbing probingGen;
4140    probingGen.setUsingObjective(1);
4141    probingGen.setMaxPass(1);
4142    probingGen.setMaxPassRoot(1);
4143    // Number of unsatisfied variables to look at
4144    probingGen.setMaxProbe(10);
4145    probingGen.setMaxProbeRoot(50);
4146    // How far to follow the consequences
4147    probingGen.setMaxLook(10);
4148    probingGen.setMaxLookRoot(50);
4149    probingGen.setMaxLookRoot(10);
4150    // Only look at rows with fewer than this number of elements
4151    probingGen.setMaxElements(200);
4152    probingGen.setMaxElementsRoot(300);
4153    probingGen.setRowCuts(3);
4154    // set default action (0=off,1=on,2=root)
4155    int probingAction=1;
4156
4157    CglKnapsackCover knapsackGen;
4158    //knapsackGen.switchOnExpensive();
4159    //knapsackGen.setMaxInKnapsack(100);
4160    // set default action (0=off,1=on,2=root)
4161    int knapsackAction=3;
4162
4163    CglRedSplit redsplitGen;
4164    //redsplitGen.setLimit(100);
4165    // set default action (0=off,1=on,2=root)
4166    // Off as seems to give some bad cuts
4167    int redsplitAction=0;
4168
4169    CglFakeClique cliqueGen(NULL,false);
4170    //CglClique cliqueGen(false,true);
4171    cliqueGen.setStarCliqueReport(false);
4172    cliqueGen.setRowCliqueReport(false);
4173    cliqueGen.setMinViolation(0.1);
4174    // set default action (0=off,1=on,2=root)
4175    int cliqueAction=3;
4176
4177    CglMixedIntegerRounding2 mixedGen;
4178    // set default action (0=off,1=on,2=root)
4179    int mixedAction=3;
4180
4181    CglFlowCover flowGen;
4182    // set default action (0=off,1=on,2=root)
4183    int flowAction=3;
4184
4185    CglTwomir twomirGen;
4186    twomirGen.setMaxElements(250);
4187    // set default action (0=off,1=on,2=root)
4188    int twomirAction=2;
4189#ifndef DEBUG_MALLOC
4190    CglLandP landpGen;
4191#endif
4192    // set default action (0=off,1=on,2=root)
4193    int landpAction=0;
4194    CglResidualCapacity residualCapacityGen;
4195    // set default action (0=off,1=on,2=root)
4196    int residualCapacityAction=0;
4197
4198#ifdef ZERO_HALF_CUTS
4199    CglZeroHalf zerohalfGen;
4200    //zerohalfGen.switchOnExpensive();
4201#endif
4202    // set default action (0=off,1=on,2=root)
4203    int zerohalfAction=0;
4204
4205    // Stored cuts
4206    bool storedCuts = false;
4207
4208    int useCosts=0;
4209    // don't use input solution
4210    int useSolution=-1;
4211   
4212    // total number of commands read
4213    int numberGoodCommands=0;
4214    // Set false if user does anything advanced
4215    bool defaultSettings=true;
4216
4217    // Hidden stuff for barrier
4218    int choleskyType = 0;
4219    int gamma=0;
4220    int scaleBarrier=0;
4221    int doKKT=0;
4222    int crossover=2; // do crossover unless quadratic
4223    // For names
4224    int lengthName = 0;
4225    std::vector<std::string> rowNames;
4226    std::vector<std::string> columnNames;
4227    // Default strategy stuff
4228    {
4229      // try changing tolerance at root
4230#define MORE_CUTS
4231#ifdef MORE_CUTS
4232      gomoryGen.setAwayAtRoot(0.005);
4233      twomirGen.setAwayAtRoot(0.005);
4234      twomirGen.setAway(0.01);
4235#else
4236      gomoryGen.setAwayAtRoot(0.01);
4237      twomirGen.setAwayAtRoot(0.01);
4238      twomirGen.setAway(0.01);
4239#endif
4240      int iParam;
4241      iParam = whichParam(DIVEOPT,numberParameters_,parameters_);
4242      parameters_[iParam].setIntValue(3);
4243      iParam = whichParam(FPUMPITS,numberParameters_,parameters_);
4244      parameters_[iParam].setIntValue(30);
4245      iParam = whichParam(FPUMPTUNE,numberParameters_,parameters_);
4246      parameters_[iParam].setIntValue(1005043);
4247      iParam = whichParam(PROCESSTUNE,numberParameters_,parameters_);
4248      parameters_[iParam].setIntValue(6);
4249      tunePreProcess=6;
4250      iParam = whichParam(DIVINGC,numberParameters_,parameters_);
4251      parameters_[iParam].setCurrentOption("on");
4252      iParam = whichParam(RINS,numberParameters_,parameters_);
4253      parameters_[iParam].setCurrentOption("on");
4254      iParam = whichParam(PROBINGCUTS,numberParameters_,parameters_);
4255      parameters_[iParam].setCurrentOption("on");
4256      probingAction=1;
4257    }
4258    std::string field;
4259    if (!noPrinting_) {
4260      sprintf(generalPrint,"Coin Cbc and Clp Solver version %s, build %s",
4261              CBCVERSION,__DATE__);
4262      generalMessageHandler->message(CLP_GENERAL,generalMessages)
4263        << generalPrint
4264        <<CoinMessageEol;
4265      // Print command line
4266      if (argc>1) {
4267        bool foundStrategy=false;
4268        sprintf(generalPrint,"command line - ");
4269        for (int i=0;i<argc;i++) {
4270          if (!argv[i])
4271            break;
4272          if (strstr(argv[i],"strat"))
4273            foundStrategy=true;
4274          sprintf(generalPrint+strlen(generalPrint),"%s ",argv[i]);
4275        }
4276        if (!foundStrategy)
4277          sprintf(generalPrint+strlen(generalPrint),"(default strategy 1)");
4278        generalMessageHandler->message(CLP_GENERAL,generalMessages)
4279          << generalPrint
4280          <<CoinMessageEol;
4281      }
4282    }
4283    while (1) {
4284      // next command
4285      field=CoinReadGetCommand(argc,argv);
4286      // Reset time
4287      time1 = CoinCpuTime();
4288      // adjust field if has odd trailing characters
4289      char temp [200];
4290      strcpy(temp,field.c_str());
4291      int length = strlen(temp);
4292      for (int k=length-1;k>=0;k--) {
4293        if (temp[k]<' ')
4294          length--;
4295        else
4296          break;
4297      }
4298      temp[length]='\0';
4299      field=temp;
4300      // exit if null or similar
4301      if (!field.length()) {
4302        if (numberGoodCommands==1&&goodModel) {
4303          // we just had file name - do branch and bound
4304          field="branch";
4305        } else if (!numberGoodCommands) {
4306          // let's give the sucker a hint
4307          std::cout
4308            <<"CoinSolver takes input from arguments ( - switches to stdin)"
4309            <<std::endl
4310            <<"Enter ? for list of commands or help"<<std::endl;
4311          field="-";
4312        } else {
4313          break;
4314        }
4315      }
4316     
4317      // see if ? at end
4318      int numberQuery=0;
4319      if (field!="?"&&field!="???") {
4320        int length = field.length();
4321        int i;
4322        for (i=length-1;i>0;i--) {
4323          if (field[i]=='?') 
4324            numberQuery++;
4325          else
4326            break;
4327        }
4328        field=field.substr(0,length-numberQuery);
4329      }
4330      // find out if valid command
4331      int iParam;
4332      int numberMatches=0;
4333      int firstMatch=-1;
4334      for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4335        int match = parameters_[iParam].matches(field);
4336        if (match==1) {
4337          numberMatches = 1;
4338          firstMatch=iParam;
4339          break;
4340        } else {
4341          if (match&&firstMatch<0)
4342            firstMatch=iParam;
4343          numberMatches += match>>1;
4344        }
4345      }
4346      if (iParam<numberParameters_&&!numberQuery) {
4347        // found
4348        CbcOrClpParam found = parameters_[iParam];
4349        CbcOrClpParameterType type = found.type();
4350        int valid;
4351        numberGoodCommands++;
4352        if (type==BAB&&goodModel) {
4353          // check if any integers
4354#ifndef CBC_OTHER_SOLVER
4355#ifdef COIN_HAS_ASL
4356          if (info.numberSos&&doSOS&&statusUserFunction_[0]) {
4357            // SOS
4358            numberSOS = info.numberSos;
4359          }
4360#endif
4361          lpSolver = clpSolver->getModelPtr();
4362          if (!lpSolver->integerInformation()&&!numberSOS&&
4363              !clpSolver->numberSOS()&&!model_.numberObjects()&&!clpSolver->numberObjects())
4364            type=DUALSIMPLEX;
4365#endif
4366        }
4367        if (type==GENERALQUERY) {
4368          bool evenHidden=false;
4369          if ((verbose&8)!=0) {
4370            // even hidden
4371            evenHidden = true;
4372            verbose &= ~8;
4373          }
4374#ifdef COIN_HAS_ASL
4375          if (verbose<4&&statusUserFunction_[0])
4376            verbose +=4;
4377#endif
4378          if (verbose<4) {
4379            std::cout<<"In argument list keywords have leading - "
4380              ", -stdin or just - switches to stdin"<<std::endl;
4381            std::cout<<"One command per line (and no -)"<<std::endl;
4382            std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
4383            std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
4384            std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
4385            std::cout<<"abcd value sets value"<<std::endl;
4386            std::cout<<"Commands are:"<<std::endl;
4387          } else {
4388            std::cout<<"Cbc options are set within AMPL with commands like:"<<std::endl<<std::endl;
4389            std::cout<<"         option cbc_options \"cuts=root log=2 feas=on slog=1\""<<std::endl<<std::endl;
4390            std::cout<<"only maximize, dual, primal, help and quit are recognized without ="<<std::endl;
4391          }
4392          int maxAcross=5;
4393          if ((verbose%4)!=0)
4394            maxAcross=1;
4395          int limits[]={1,51,101,151,201,251,301,351,401};
4396          std::vector<std::string> types;
4397          types.push_back("Double parameters:");
4398          types.push_back("Branch and Cut double parameters:");
4399          types.push_back("Integer parameters:");
4400          types.push_back("Branch and Cut integer parameters:");
4401          types.push_back("Keyword parameters:");
4402          types.push_back("Branch and Cut keyword parameters:");
4403          types.push_back("Actions or string parameters:");
4404          types.push_back("Branch and Cut actions:");
4405          int iType;
4406          for (iType=0;iType<8;iType++) {
4407            int across=0;
4408            if ((verbose%4)!=0)
4409              std::cout<<std::endl;
4410            std::cout<<types[iType]<<std::endl;
4411            if ((verbose&2)!=0)
4412              std::cout<<std::endl;
4413            for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4414              int type = parameters_[iParam].type();
4415              if ((parameters_[iParam].displayThis()||evenHidden)&&
4416                  type>=limits[iType]
4417                  &&type<limits[iType+1]) {
4418                // but skip if not useful for ampl (and in ampl mode)
4419                if (verbose>=4&&(parameters_[iParam].whereUsed()&4)==0)
4420                  continue;
4421                if (!across) {
4422                  if ((verbose&2)==0) 
4423                    std::cout<<"  ";
4424                  else
4425                    std::cout<<"Command ";
4426                }
4427                std::cout<<parameters_[iParam].matchName()<<"  ";
4428                across++;
4429                if (across==maxAcross) {
4430                  across=0;
4431                  if ((verbose%4)!=0) {
4432                    // put out description as well
4433                    if ((verbose&1)!=0) 
4434                      std::cout<<parameters_[iParam].shortHelp();
4435                    std::cout<<std::endl;
4436                    if ((verbose&2)!=0) {
4437                      std::cout<<"---- description"<<std::endl;
4438                      parameters_[iParam].printLongHelp();
4439                      std::cout<<"----"<<std::endl<<std::endl;
4440                    }
4441                  } else {
4442                    std::cout<<std::endl;
4443                  }
4444                }
4445              }
4446            }
4447            if (across)
4448              std::cout<<std::endl;
4449          }
4450        } else if (type==FULLGENERALQUERY) {
4451          std::cout<<"Full list of commands is:"<<std::endl;
4452          int maxAcross=5;
4453          int limits[]={1,51,101,151,201,251,301,351,401};
4454          std::vector<std::string> types;
4455          types.push_back("Double parameters:");
4456          types.push_back("Branch and Cut double parameters:");
4457          types.push_back("Integer parameters:");
4458          types.push_back("Branch and Cut integer parameters:");
4459          types.push_back("Keyword parameters:");
4460          types.push_back("Branch and Cut keyword parameters:");
4461          types.push_back("Actions or string parameters:");
4462          types.push_back("Branch and Cut actions:");
4463          int iType;
4464          for (iType=0;iType<8;iType++) {
4465            int across=0;
4466            std::cout<<types[iType]<<"  ";
4467            for ( iParam=0; iParam<numberParameters_; iParam++ ) {
4468              int type = parameters_[iParam].type();
4469              if (type>=limits[iType]
4470                  &&type<limits[iType+1]) {
4471                if (!across)
4472                  std::cout<<"  ";
4473                std::cout<<parameters_[iParam].matchName()<<"  ";
4474                across++;
4475                if (across==maxAcross) {
4476                  std::cout<<std::endl;
4477                  across=0;
4478                }
4479              }
4480            }
4481            if (across)
4482              std::cout<<std::endl;
4483          }
4484        } else if (type<101) {
4485          // get next field as double
4486          double value = CoinReadGetDoubleField(argc,argv,&valid);
4487          if (!valid) {
4488            if (type<51) {
4489              int returnCode;
4490              const char * message = 
4491                parameters_[iParam].setDoubleParameterWithMessage(lpSolver,value,returnCode);
4492              if (!noPrinting_&&strlen(message)) {
4493                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4494                  << message
4495                  <<CoinMessageEol;
4496              }
4497            } else if (type<81) {
4498              int returnCode;
4499              const char * message = 
4500                parameters_[iParam].setDoubleParameterWithMessage(model_,value,returnCode);
4501              if (!noPrinting_&&strlen(message)) {
4502                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4503                  << message
4504                  <<CoinMessageEol;
4505              }
4506            } else {
4507              int returnCode;
4508              const char * message = 
4509                parameters_[iParam].setDoubleParameterWithMessage(lpSolver,value,returnCode);
4510              if (!noPrinting_&&strlen(message)) {
4511                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4512                  << message
4513                  <<CoinMessageEol;
4514              }
4515              switch(type) {
4516              case DJFIX:
4517                djFix=value;
4518#ifndef CBC_OTHER_SOLVER
4519                if (goodModel&&djFix<1.0e20) {
4520                  // do some fixing
4521                  clpSolver = dynamic_cast< OsiClpSolverInterface*> (model_.solver());
4522                  clpSolver->initialSolve();
4523                  lpSolver = clpSolver->getModelPtr();
4524                  int numberColumns = lpSolver->numberColumns();
4525                  int i;
4526                  const char * type = lpSolver->integerInformation();
4527                  double * lower = lpSolver->columnLower();
4528                  double * upper = lpSolver->columnUpper();
4529                  double * solution = lpSolver->primalColumnSolution();
4530                  double * dj = lpSolver->dualColumnSolution();
4531                  int numberFixed=0;
4532                  double dextra4 = parameters_[whichParam(DEXTRA4,numberParameters_,parameters_)].doubleValue();
4533                  if (dextra4)
4534                    printf("Multiple for continuous dj fixing is %g\n",dextra4);
4535                  for (i=0;i<numberColumns;i++) {
4536                    double djValue = dj[i];
4537                    if (!type[i])
4538                      djValue *= dextra4;
4539                    if (type[i]||dextra4) {
4540                      double value = solution[i];
4541                      if (value<lower[i]+1.0e-5&&djValue>djFix) {
4542                        solution[i]=lower[i];
4543                        upper[i]=lower[i];
4544                        numberFixed++;
4545                      } else if (value>upper[i]-1.0e-5&&djValue<-djFix) {
4546                        solution[i]=upper[i];
4547                        lower[i]=upper[i];
4548                        numberFixed++;
4549                      }
4550                    }
4551                  }
4552                  sprintf(generalPrint,"%d columns fixed\n",numberFixed);
4553                  generalMessageHandler->message(CLP_GENERAL,generalMessages)
4554                    << generalPrint
4555                    <<CoinMessageEol;
4556                }
4557#endif
4558                break;
4559              case TIGHTENFACTOR:
4560                tightenFactor=value;
4561                if(!complicatedInteger)
4562                  defaultSettings=false; // user knows what she is doing
4563                break;
4564              default:
4565                break;
4566              }
4567            }
4568          } else if (valid==1) {
4569            std::cout<<" is illegal for double parameter "<<parameters_[iParam].name()<<" value remains "<<
4570              parameters_[iParam].doubleValue()<<std::endl;
4571          } else {
4572            std::cout<<parameters_[iParam].name()<<" has value "<<
4573              parameters_[iParam].doubleValue()<<std::endl;
4574          }
4575        } else if (type<201) {
4576          // get next field as int
4577          int value = CoinReadGetIntField(argc,argv,&valid);
4578          if (!valid) {
4579            if (type<151) {
4580              if (parameters_[iParam].type()==PRESOLVEPASS)
4581                preSolve = value;
4582              else if (parameters_[iParam].type()==IDIOT)
4583                doIdiot = value;
4584              else if (parameters_[iParam].type()==SPRINT)
4585                doSprint = value;
4586              else if (parameters_[iParam].type()==OUTPUTFORMAT)
4587                outputFormat = value;
4588              else if (parameters_[iParam].type()==SLPVALUE)
4589                slpValue = value;
4590              else if (parameters_[iParam].type()==CPP)
4591                cppValue = value;
4592              else if (parameters_[iParam].type()==PRESOLVEOPTIONS)
4593                presolveOptions = value;
4594              else if (parameters_[iParam].type()==PRINTOPTIONS)
4595                printOptions = value;
4596              else if (parameters_[iParam].type()==SUBSTITUTION)
4597                substitution = value;
4598              else if (parameters_[iParam].type()==DUALIZE)
4599                dualize = value;
4600              else if (parameters_[iParam].type()==PROCESSTUNE)
4601                tunePreProcess = value;
4602              else if (parameters_[iParam].type()==USESOLUTION)
4603                useSolution = value;
4604              else if (parameters_[iParam].type()==VERBOSE)
4605                verbose = value;
4606              int returnCode;
4607              const char * message = 
4608                parameters_[iParam].setIntParameterWithMessage(lpSolver,value,returnCode);
4609              if (!noPrinting_&&strlen(message)) {
4610                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4611                  << message
4612                  <<CoinMessageEol;
4613              }
4614            } else {
4615              if (parameters_[iParam].type()==CUTPASS)
4616                cutPass = value;
4617              else if (parameters_[iParam].type()==CUTPASSINTREE)
4618                cutPassInTree = value;
4619              else if (parameters_[iParam].type()==STRONGBRANCHING||
4620                       parameters_[iParam].type()==NUMBERBEFORE)
4621                strongChanged=true;
4622              else if (parameters_[iParam].type()==FPUMPTUNE||
4623                       parameters_[iParam].type()==FPUMPITS)
4624                pumpChanged=true;
4625              else if (parameters_[iParam].type()==EXPERIMENT) {
4626                if (value>=1) {
4627                  if (!noPrinting_) {
4628                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4629                      <<"switching on dense factorization if small, and maybe fast fathoming"
4630                      <<CoinMessageEol;
4631                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4632                      <<"Gomory cuts using tolerance of 0.01 at root"
4633                      <<CoinMessageEol;
4634                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4635                      <<"Possible restart after 100 nodes if can fix many"
4636                      <<CoinMessageEol;
4637                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
4638                      <<"extra options - -diving C -diveopt 3 -rins on -tune 6 -probing on -passf 30!"
4639                      <<CoinMessageEol;
4640                  }
4641                  // try changing tolerance at root
4642                  //gomoryGen.setAwayAtRoot(0.01);
4643                  int iParam;
4644                  iParam = whichParam(DIVEOPT,numberParameters_,parameters_);
4645                  parameters_[iParam].setIntValue(3);
4646                  iParam = whichParam(FPUMPITS,numberParameters_,parameters_);
4647                  parameters_[iParam].setIntValue(30);
4648                  iParam = whichParam(FPUMPTUNE,numberParameters_,parameters_);
4649                  parameters_[iParam].setIntValue(1005043);
4650                  iParam = whichParam(PROCESSTUNE,numberParameters_,parameters_);
4651                  parameters_[iParam].setIntValue(6);
4652                  tunePreProcess=6;
4653                  iParam = whichParam(DIVINGC,numberParameters_,parameters_);
4654                  parameters_[iParam].setCurrentOption("on");
4655                  iParam = whichParam(RINS,numberParameters_,parameters_);
4656                  parameters_[iParam].setCurrentOption("on");
4657                  iParam = whichParam(PROBINGCUTS,numberParameters_,parameters_);
4658                  parameters_[iParam].setCurrentOption("on");
4659                  probingAction=1;
4660                }
4661              } else if (parameters_[iParam].type()==STRATEGY) {
4662                if (value==0) {
4663                  gomoryGen.setAwayAtRoot(0.05);
4664                  int iParam;
4665                  iParam = whichParam(DIVEOPT,numberParameters_,parameters_);
4666                  parameters_[iParam].setIntValue(-1);
4667                  iParam = whichParam(FPUMPITS,numberParameters_,parameters_);
4668                  parameters_[iParam].setIntValue(20);
4669                  iParam = whichParam(FPUMPTUNE,numberParameters_,parameters_);
4670                  parameters_[iParam].setIntValue(1003);
4671                  iParam = whichParam(PROCESSTUNE,numberParameters_,parameters_);
4672                  parameters_[iParam].setIntValue(-1);
4673                  tunePreProcess=0;
4674                  iParam = whichParam(DIVINGC,numberParameters_,parameters_);
4675                  parameters_[iParam].setCurrentOption("off");
4676                  iParam = whichParam(RINS,numberParameters_,parameters_);
4677                  parameters_[iParam].setCurrentOption("off");
4678                  iParam = whichParam(PROBINGCUTS,numberParameters_,parameters_);
4679                  parameters_[iParam].setCurrentOption("on");
4680                  probingAction=1;
4681                }
4682              }
4683              int returnCode;
4684              const char * message = 
4685                parameters_[iParam].setIntParameterWithMessage(model_,value,returnCode);
4686              if (!noPrinting_&&strlen(message)) {
4687                generalMessageHandler->message(CLP_GENERAL,generalMessages)
4688                  << message
4689                  <<CoinMessageEol;
4690              }
4691            }
4692          } else if (valid==1) {
4693            std::cout<<" is illegal for integer parameter "<<parameters_[iParam].name()<<" value remains "<<
4694              parameters_[iParam].intValue()<<std::endl;
4695          } else {
4696            std::cout<<parameters_[iParam].name()<<" has value "<<
4697              parameters_[iParam].intValue()<<std::endl;
4698          }
4699        } else if (type<301) {
4700          // one of several strings
4701          std::string value = CoinReadGetString(argc,argv);
4702          int action = parameters_[iParam].parameterOption(value);
4703          if (action<0) {
4704            if (value!="EOL") {
4705              // no match
4706              parameters_[iParam].printOptions();
4707            } else {
4708              // print current value
4709              std::cout<<parameters_[iParam].name()<<" has value "<<
4710                parameters_[iParam].currentOption()<<std::endl;
4711            }
4712          } else {
4713            const char * message = 
4714              parameters_[iParam].setCurrentOptionWithMessage(action);
4715            if (!noPrinting_&&strlen(message)) {
4716              generalMessageHandler->message(CLP_GENERAL,generalMessages)
4717                << message
4718                <<CoinMessageEol;
4719            }
4720            // for now hard wired
4721            switch (type) {
4722            case DIRECTION:
4723              if (action==0)
4724                lpSolver->setOptimizationDirection(1);
4725              else if (action==1)
4726                lpSolver->setOptimizationDirection(-1);
4727              else
4728                lpSolver->setOptimizationDirection(0);
4729              break;
4730            case DUALPIVOT:
4731              if (action==0) {
4732                ClpDualRowSteepest steep(3);
4733                lpSolver->setDualRowPivotAlgorithm(steep);
4734              } else if (action==1) {
4735                ClpDualRowDantzig dantzig;
4736                //ClpDualRowSteepest dantzig(5);
4737                lpSolver->setDualRowPivotAlgorithm(dantzig);
4738              } else if (action==2) {
4739                // partial steep
4740                ClpDualRowSteepest steep(2);
4741                lpSolver->setDualRowPivotAlgorithm(steep);
4742              } else {
4743                ClpDualRowSteepest steep;
4744                lpSolver->setDualRowPivotAlgorithm(steep);
4745              }
4746              break;
4747            case PRIMALPIVOT:
4748              if (action==0) {
4749                ClpPrimalColumnSteepest steep(3);
4750                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4751              } else if (action==1) {
4752                ClpPrimalColumnSteepest steep(0);
4753                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4754              } else if (action==2) {
4755                ClpPrimalColumnDantzig dantzig;
4756                lpSolver->setPrimalColumnPivotAlgorithm(dantzig);
4757              } else if (action==3) {
4758                ClpPrimalColumnSteepest steep(4);
4759                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4760              } else if (action==4) {
4761                ClpPrimalColumnSteepest steep(1);
4762                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4763              } else if (action==5) {
4764                ClpPrimalColumnSteepest steep(2);
4765                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4766              } else if (action==6) {
4767                ClpPrimalColumnSteepest steep(10);
4768                lpSolver->setPrimalColumnPivotAlgorithm(steep);
4769              }
4770              break;
4771            case SCALING:
4772              lpSolver->scaling(action);
4773              solver->setHintParam(OsiDoScale,action!=0,OsiHintTry);
4774              doScaling = action;
4775              break;
4776            case AUTOSCALE:
4777              lpSolver->setAutomaticScaling(action!=0);
4778              break;
4779            case SPARSEFACTOR:
4780              lpSolver->setSparseFactorization((1-action)!=0);
4781              break;
4782            case BIASLU:
4783              lpSolver->factorization()->setBiasLU(action);
4784              break;
4785            case PERTURBATION:
4786              if (action==0)
4787                lpSolver->setPerturbation(50);
4788              else
4789                lpSolver->setPerturbation(100);
4790              break;
4791            case ERRORSALLOWED:
4792              allowImportErrors = action;
4793              break;
4794            case INTPRINT:
4795              printMode=action;
4796              break;
4797              //case ALGORITHM:
4798              //algorithm  = action;
4799              //defaultSettings=false; // user knows what she is doing
4800              //abort();
4801              //break;
4802            case KEEPNAMES:
4803              keepImportNames = 1-action;
4804              break;
4805            case PRESOLVE:
4806              if (action==0)
4807                preSolve = 5;
4808              else if (action==1)
4809                preSolve=0;
4810              else if (action==2)
4811                preSolve=10;
4812              else
4813                preSolveFile=true;
4814              break;
4815            case PFI:
4816              lpSolver->factorization()->setForrestTomlin(action==0);
4817              break;
4818            case CRASH:
4819              doCrash=action;
4820              break;
4821            case VECTOR:
4822              doVector=action;
4823              break;
4824            case MESSAGES:
4825              lpSolver->messageHandler()->setPrefix(action!=0);
4826              break;
4827            case CHOLESKY:
4828              choleskyType = action;
4829              break;
4830            case GAMMA:
4831              gamma=action;
4832              break;
4833            case BARRIERSCALE:
4834              scaleBarrier=action;
4835              break;
4836            case KKT:
4837              doKKT=action;
4838              break;
4839            case CROSSOVER:
4840              crossover=action;
4841              break;
4842            case SOS:
4843              doSOS=action;
4844              break;
4845            case GOMORYCUTS:
4846              defaultSettings=false; // user knows what she is doing
4847              gomoryAction = action;
4848              break;
4849            case PROBINGCUTS:
4850              defaultSettings=false; // user knows what she is doing
4851              probingAction = action;
4852              break;
4853            case KNAPSACKCUTS:
4854              defaultSettings=false; // user knows what she is doing
4855              knapsackAction = action;
4856              break;
4857            case REDSPLITCUTS:
4858              defaultSettings=false; // user knows what she is doing
4859              redsplitAction = action;
4860              break;
4861            case CLIQUECUTS:
4862              defaultSettings=false; // user knows what she is doing
4863              cliqueAction = action;
4864              break;
4865            case FLOWCUTS:
4866              defaultSettings=false; // user knows what she is doing
4867              flowAction = action;
4868              break;
4869            case MIXEDCUTS:
4870              defaultSettings=false; // user knows what she is doing
4871              mixedAction = action;
4872              break;
4873            case TWOMIRCUTS:
4874              defaultSettings=false; // user knows what she is doing
4875              twomirAction = action;
4876              break;
4877            case LANDPCUTS:
4878              defaultSettings=false; // user knows what she is doing
4879              landpAction = action;
4880              break;
4881            case RESIDCUTS:
4882              defaultSettings=false; // user knows what she is doing
4883              residualCapacityAction = action;
4884              break;
4885            case ZEROHALFCUTS:
4886              defaultSettings=false; // user knows what she is doing
4887              zerohalfAction = action;
4888              break;
4889            case ROUNDING:
4890              defaultSettings=false; // user knows what she is doing
4891              break;
4892            case FPUMP:
4893              defaultSettings=false; // user knows what she is doing
4894              break;
4895            case RINS:
4896              break;
4897            case DINS:
4898              break;
4899            case RENS:
4900              break;
4901            case CUTSSTRATEGY:
4902              gomoryAction = action;
4903              probingAction = action;
4904              knapsackAction = action;
4905              zerohalfAction = action;
4906              cliqueAction = action;
4907              flowAction = action;
4908              mixedAction = action;
4909              twomirAction = action;
4910              //landpAction = action;
4911              parameters_[whichParam(GOMORYCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4912              parameters_[whichParam(PROBINGCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4913              parameters_[whichParam(KNAPSACKCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4914              parameters_[whichParam(CLIQUECUTS,numberParameters_,parameters_)].setCurrentOption(action);
4915              parameters_[whichParam(FLOWCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4916              parameters_[whichParam(MIXEDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4917              parameters_[whichParam(TWOMIRCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4918              parameters_[whichParam(ZEROHALFCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4919              if (!action) {
4920                redsplitAction = action;
4921                parameters_[whichParam(REDSPLITCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4922                landpAction = action;
4923                parameters_[whichParam(LANDPCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4924                residualCapacityAction = action;
4925                parameters_[whichParam(RESIDCUTS,numberParameters_,parameters_)].setCurrentOption(action);
4926              }
4927              break;
4928            case HEURISTICSTRATEGY:
4929              parameters_[whichParam(ROUNDING,numberParameters_,parameters_)].setCurrentOption(action);
4930              parameters_[whichParam(GREEDY,numberParameters_,parameters_)].setCurrentOption(action);
4931              parameters_[whichParam(COMBINE,numberParameters_,parameters_)].setCurrentOption(action);
4932              //parameters_[whichParam(LOCALTREE,numberParameters_,parameters_)].setCurrentOption(action);
4933              parameters_[whichParam(FPUMP,numberParameters_,parameters_)].setCurrentOption(action);
4934              parameters_[whichParam(DIVINGC,numberParameters_,parameters_)].setCurrentOption(action);
4935              parameters_[whichParam(RINS,numberParameters_,parameters_)].setCurrentOption(action);
4936              break;
4937            case GREEDY:
4938            case DIVINGS:
4939            case DIVINGC:
4940            case DIVINGF:
4941            case DIVINGG:
4942            case DIVINGL:
4943            case DIVINGP:
4944            case DIVINGV:
4945            case COMBINE:
4946            case PIVOTANDFIX:
4947            case RANDROUND:
4948            case LOCALTREE:
4949            case NAIVE:
4950            case CPX:
4951              defaultSettings=false; // user knows what she is doing
4952              break;
4953            case COSTSTRATEGY:
4954              useCosts=action;
4955              break;
4956            case NODESTRATEGY:
4957              nodeStrategy=action;
4958              break;
4959            case PREPROCESS:
4960              preProcess = action;
4961              break;
4962              break;
4963            default:
4964              abort();
4965            }
4966          }
4967        } else {
4968          // action
4969          if (type==EXIT) {
4970#ifdef COIN_HAS_ASL
4971            if(statusUserFunction_[0]) {
4972              if (info.numberIntegers||info.numberBinary) {
4973                // integer
4974              } else {
4975                // linear
4976              }
4977              writeAmpl(&info);
4978              freeArrays2(&info);
4979              freeArgs(&info);
4980            }
4981#endif
4982            break; // stop all
4983          }
4984          switch (type) {
4985          case DUALSIMPLEX:
4986          case PRIMALSIMPLEX:
4987          case SOLVECONTINUOUS:
4988          case BARRIER:
4989            if (goodModel) {
4990              double objScale = 
4991                parameters_[whichParam(OBJSCALE2,numberParameters_,parameters_)].doubleValue();
4992              if (objScale!=1.0) {
4993                int iColumn;
4994                int numberColumns=lpSolver->numberColumns();
4995                double * dualColumnSolution = 
4996                  lpSolver->dualColumnSolution();
4997                ClpObjective * obj = lpSolver->objectiveAsObject();
4998                assert(dynamic_cast<ClpLinearObjective *> (obj));
4999                double offset;
5000                double * objective = obj->gradient(NULL,NULL,offset,true);
5001                for (iColumn=0;iColumn<numberColumns;iColumn++) {
5002                  dualColumnSolution[iColumn] *= objScale;
5003                  objective[iColumn] *= objScale;;
5004                }
5005                int iRow;
5006                int numberRows=lpSolver->numberRows();
5007                double * dualRowSolution = 
5008                  lpSolver->dualRowSolution();
5009                for (iRow=0;iRow<numberRows;iRow++) 
5010                  dualRowSolution[iRow] *= objScale;
5011                lpSolver->setObjectiveOffset(objScale*lpSolver->objectiveOffset());
5012              }
5013              ClpSolve::SolveType method;
5014              ClpSolve::PresolveType presolveType;
5015              ClpSimplex * model2 = lpSolver;
5016              if (dualize) {
5017                bool tryIt=true;
5018                double fractionColumn=1.0;
5019                double fractionRow=1.0;
5020                if (dualize==3) {
5021                  dualize=1;
5022                  int numberColumns=lpSolver->numberColumns();
5023                  int numberRows=lpSolver->numberRows();
5024                  if (numberRows<50000||5*numberColumns>numberRows) {
5025                    tryIt=false;
5026                  } else {
5027                    fractionColumn=0.1;
5028                    fractionRow=0.1;
5029                  }
5030                }
5031                if (tryIt) {
5032                  model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(fractionRow,fractionColumn);
5033                  if (model2) {
5034                    sprintf(generalPrint,"Dual of model has %d rows and %d columns",
5035                            model2->numberRows(),model2->numberColumns());
5036                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
5037                      << generalPrint
5038                      <<CoinMessageEol;
5039                    model2->setOptimizationDirection(1.0);
5040                  } else {
5041                    model2 = lpSolver;
5042                    dualize=0;
5043                  }
5044                } else {
5045                  dualize=0;
5046                }
5047              }
5048              if (noPrinting_)
5049                lpSolver->setLogLevel(0);
5050              ClpSolve solveOptions;
5051              solveOptions.setPresolveActions(presolveOptions);
5052              solveOptions.setSubstitution(substitution);
5053              if (preSolve!=5&&preSolve) {
5054                presolveType=ClpSolve::presolveNumber;
5055                if (preSolve<0) {
5056                  preSolve = - preSolve;
5057                  if (preSolve<=100) {
5058                    presolveType=ClpSolve::presolveNumber;
5059                    sprintf(generalPrint,"Doing %d presolve passes - picking up non-costed slacks",
5060                           preSolve);
5061                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
5062                      << generalPrint
5063                      <<CoinMessageEol;
5064                    solveOptions.setDoSingletonColumn(true);
5065                  } else {
5066                    preSolve -=100;
5067                    presolveType=ClpSolve::presolveNumberCost;
5068                    sprintf(generalPrint,"Doing %d presolve passes - picking up costed slacks",
5069                           preSolve);
5070                    generalMessageHandler->message(CLP_GENERAL,generalMessages)
5071                      << generalPrint
5072                      <<CoinMessageEol;
5073                  }
5074                } 
5075              } else if (preSolve) {
5076                presolveType=ClpSolve::presolveOn;
5077              } else {
5078                presolveType=ClpSolve::presolveOff;
5079              }
5080              solveOptions.setPresolveType(presolveType,preSolve);
5081              if (type==DUALSIMPLEX||type==SOLVECONTINUOUS) {
5082                method=ClpSolve::useDual;
5083              } else if (type==PRIMALSIMPLEX) {
5084                method=ClpSolve::usePrimalorSprint;
5085              } else {
5086                method = ClpSolve::useBarrier;
5087                if (crossover==1) {
5088                  method=ClpSolve::useBarrierNoCross;
5089                } else if (crossover==2) {
5090                  ClpObjective * obj = lpSolver->objectiveAsObject();
5091                  if (obj->type()>1) {
5092                    method=ClpSolve::useBarrierNoCross;
5093                    presolveType=ClpSolve::presolveOff;
5094                    solveOptions.setPresolveType(presolveType,preSolve);
5095                  } 
5096                }
5097              }
5098              solveOptions.setSolveType(method);
5099              if(preSolveFile)
5100                presolveOptions |= 0x40000000;
5101              solveOptions.setSpecialOption(4,presolveOptions);
5102              solveOptions.setSpecialOption(5,printOptions);
5103              if (doVector) {
5104                ClpMatrixBase * matrix = lpSolver->clpMatrix();
5105                if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
5106                  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
5107                  clpMatrix->makeSpecialColumnCopy();
5108                }
5109              }
5110              if (method==ClpSolve::useDual) {
5111                // dual
5112                if (doCrash)
5113                  solveOptions.setSpecialOption(0,1,doCrash); // crash
5114                else if (doIdiot)
5115                  solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
5116              } else if (method==ClpSolve::usePrimalorSprint) {
5117                // primal
5118                // if slp turn everything off
5119                if (slpValue>0) {
5120                  doCrash=false;
5121                  doSprint=0;
5122                  doIdiot=-1;
5123                  solveOptions.setSpecialOption(1,10,slpValue); // slp
5124                  method=ClpSolve::usePrimal;
5125                }
5126                if (doCrash) {
5127                  solveOptions.setSpecialOption(1,1,doCrash); // crash
5128                } else if (doSprint>0) {
5129                  // sprint overrides idiot
5130                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
5131                } else if (doIdiot>0) {
5132                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
5133                } else if (slpValue<=0) {
5134                  if (doIdiot==0) {
5135                    if (doSprint==0)
5136                      solveOptions.setSpecialOption(1,4); // all slack
5137                    else
5138                      solveOptions.setSpecialOption(1,9); // all slack or sprint
5139                  } else {
5140                    if (doSprint==0)
5141                      solveOptions.setSpecialOption(1,8); // all slack or idiot
5142                    else
5143                      solveOptions.setSpecialOption(1,7); // initiative
5144                  }
5145                }
5146                if (basisHasValues==-1)
5147                  solveOptions.setSpecialOption(1,11); // switch off values
5148              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
5149                int barrierOptions = choleskyType;
5150                if (scaleBarrier)
5151                  barrierOptions |= 8;
5152                if (doKKT)
5153                  barrierOptions |= 16;
5154                if (gamma)
5155                  barrierOptions |= 32*gamma;
5156                if (crossover==3) 
5157                  barrierOptions |= 256; // try presolve in crossover
5158                solveOptions.setSpecialOption(4,barrierOptions);
5159              }
5160              model2->setMaximumSeconds(model_.getMaximumSeconds());
5161#ifdef COIN_HAS_LINK
5162              OsiSolverInterface * coinSolver = model_.solver();
5163              OsiSolverLink * linkSolver = dynamic_cast< OsiSolverLink*> (coinSolver);
5164              if (!linkSolver) {
5165                model2->initialSolve(solveOptions);
5166              } else {
5167                // special solver
5168                int testOsiOptions = parameters_[whichParam(TESTOSI,numberParameters_,parameters_)].intValue();
5169                double * solution = NULL;
5170                if (testOsiOptions<10) {
5171                  solution = linkSolver->nonlinearSLP(slpValue>0 ? slpValue : 20 ,1.0e-5);
5172                } else if (testOsiOptions>=10) {
5173                  CoinModel coinModel = *linkSolver->coinModel();
5174                  ClpSimplex * tempModel = approximateSolution(coinModel,slpValue>0 ? slpValue : 50 ,1.0e-5,0);
5175                  assert (tempModel);
5176                  solution = CoinCopyOfArray(tempModel->primalColumnSolution(),coinModel.numberColumns());
5177                  model2->setObjectiveValue(tempModel->objectiveValue());
5178                  model2->setProblemStatus(tempModel->problemStatus());
5179                  model2->setSecondaryStatus(tempModel->secondaryStatus());
5180                  delete tempModel;
5181                }
5182                if (solution) {
5183                  memcpy(model2->primalColumnSolution(),solution,
5184                         CoinMin(model2->numberColumns(),linkSolver->coinModel()->numberColumns())*sizeof(double));
5185                  delete [] solution;
5186                } else {
5187                  printf("No nonlinear solution\n");
5188                }
5189              }
5190#else
5191              model2->initialSolve(solveOptions);
5192#endif
5193              {
5194                // map states
5195                /* clp status
5196                   -1 - unknown e.g. before solve or if postSolve says not optimal
5197                   0 - optimal
5198                   1 - primal infeasible
5199                   2 - dual infeasible
5200                   3 - stopped on iterations or time
5201                   4 - stopped due to errors
5202                   5 - stopped by event handler (virtual int ClpEventHandler::event()) */
5203                /* cbc status
5204                   -1 before branchAndBound
5205                   0 finished - check isProvenOptimal or isProvenInfeasible to see if solution found
5206                   (or check value of best solution)
5207                   1 stopped - on maxnodes, maxsols, maxtime
5208                   2 difficulties so run was abandoned
5209                   (5 event user programmed event occurred) */
5210                /* clp secondary status of problem - may get extended
5211                   0 - none
5212                   1 - primal infeasible because dual limit reached OR probably primal
5213                   infeasible but can't prove it (main status 4)
5214                   2 - scaled problem optimal - unscaled problem has primal infeasibilities
5215                   3 - scaled problem optimal - unscaled problem has dual infeasibilities
5216                   4 - scaled problem optimal - unscaled problem has primal and dual infeasibilities
5217                   5 - giving up in primal with flagged variables
5218                   6 - failed due to empty problem check
5219                   7 - postSolve says not optimal
5220                   8 - failed due to bad element check
5221                   9 - status was 3 and stopped on time
5222                   100 up - translation of enum from ClpEventHandler
5223                */
5224                /* cbc secondary status of problem
5225                   -1 unset (status_ will also be -1)
5226                   0 search completed with solution
5227                   1 linear relaxation not feasible (or worse than cutoff)
5228                   2 stopped on gap
5229                   3 stopped on nodes
5230                   4 stopped on time
5231                   5 stopped on user event
5232                   6 stopped on solutions
5233                   7 linear relaxation unbounded
5234                */
5235                int iStatus = model2->status();
5236                int iStatus2 = model2->secondaryStatus();
5237                if (iStatus==0) {
5238                  iStatus2=0;
5239                } else if (iStatus==1) {
5240                  iStatus=0;
5241                  iStatus2=1; // say infeasible
5242                } else if (iStatus==2) {
5243                  iStatus=0;
5244                  iStatus2=7; // say unbounded
5245                } else if (iStatus==3) {
5246                  iStatus=1;
5247                  if (iStatus2==9)
5248                    iStatus2=4;
5249                  else
5250                    iStatus2=3; // Use nodes - as closer than solutions
5251                } else if (iStatus==4) {
5252                  iStatus=2; // difficulties
5253                  iStatus2=0; 
5254                }
5255                model_.setProblemStatus(iStatus);
5256                model_.setSecondaryStatus(iStatus2);
5257                //assert (lpSolver==clpSolver->getModelPtr());
5258                assert (clpSolver==model_.solver());
5259                clpSolver->setWarmStart(NULL);
5260                // and in babModel if exists
5261                if (babModel_) {
5262                  babModel_->setProblemStatus(iStatus);
5263                  babModel_->setSecondaryStatus(iStatus2);
5264                } 
5265#if NEW_STYLE_SOLVER
5266                int returnCode = callBack_->callBack(&model_,1);
5267#else
5268                int returnCode=callBack(&model,1);
5269#endif
5270                if (returnCode) {
5271                  // exit if user wants
5272#if NEW_STYLE_SOLVER
5273                  updateModel(model2,returnMode);
5274#else
5275                  delete babModel_;
5276                  babModel_ = NULL;
5277#endif
5278                  return returnCode;
5279                }
5280              }
5281              basisHasValues=1;
5282              if (dualize) {
5283                int returnCode=static_cast<ClpSimplexOther *> (lpSolver)->restoreFromDual(model2);
5284                if (model2->status()==3)
5285                  returnCode=0;
5286                delete model2;
5287                if (returnCode&&dualize!=2)
5288                  lpSolver->primal(1);
5289                model2=lpSolver;
5290              }
5291#if NEW_STYLE_SOLVER
5292              updateModel(model2,returnMode);
5293              for (iUser=0;iUser<numberUserFunctions_;iUser++) {
5294                if (statusUserFunction_[iUser])
5295                  userFunction_[iUser]->exportSolution(this,1);
5296              }
5297#endif
5298#ifdef COIN_HAS_ASL
5299              if (statusUserFunction_[0]) {
5300                double value = model2->getObjValue()*model2->getObjSense();
5301                char buf[300];
5302                int pos=0;
5303                int iStat = model2->status();
5304                if (iStat==0) {
5305                  pos +=</