source: branches/sandbox/Cbc/src/CbcSolver.cpp @ 1373

Last change on this file since 1373 was 1373, checked in by bjarni, 10 years ago

Renamed parameter constants in CbcParam?, CbcSolver?, and ClpAmplStuff? with same names as CbcOrClpParam? in CLP to make them more readable/search-able

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