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

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

changes to try and improve performance

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