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

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

final changes before cleaning

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