source: branches/reengAnn/Cbc/src/CbcSolver.cpp @ 1304

Last change on this file since 1304 was 1304, checked in by lou, 10 years ago

Merge Cbc/trunk r1262:1270 to sync with stable/2.4 (created at r1270).

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