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

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

Removed all code related to NEW_STYLE_SOLVER from CbcSolver?. A few minor
comments added to CbcModel?.

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