source: trunk/Cbc/src/CbcSolverHeuristics.cpp @ 1573

Last change on this file since 1573 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

File size: 73.5 KB
Line 
1/* $Id$ */
2// Copyright (C) 2007, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7/*! \file CbcSolverHeuristics.cpp
8    \brief Second level routines for the cbc stand-alone solver.
9*/
10
11#include "CbcConfig.h"
12#include "CoinPragma.hpp"
13
14
15#include "CoinTime.hpp"
16
17#include "OsiClpSolverInterface.hpp"
18
19#include "ClpPresolve.hpp"
20
21#include "CbcOrClpParam.hpp"
22
23#include "CbcModel.hpp"
24
25#include "CbcHeuristicLocal.hpp"
26#include "CbcHeuristicPivotAndFix.hpp"
27//#include "CbcHeuristicPivotAndComplement.hpp"
28#include "CbcHeuristicRandRound.hpp"
29#include "CbcHeuristicGreedy.hpp"
30#include "CbcHeuristicFPump.hpp"
31#include "CbcHeuristicRINS.hpp"
32
33#include "CbcHeuristicDiveCoefficient.hpp"
34#include "CbcHeuristicDiveFractional.hpp"
35#include "CbcHeuristicDiveGuided.hpp"
36#include "CbcHeuristicDiveVectorLength.hpp"
37#include "CbcHeuristicDivePseudoCost.hpp"
38#include "CbcHeuristicDiveLineSearch.hpp"
39
40#include "CbcStrategy.hpp"
41#include "OsiAuxInfo.hpp"
42
43#include "ClpSimplexOther.hpp"
44
45// Crunch down model
46void
47crunchIt(ClpSimplex * model)
48{
49#ifdef JJF_ZERO
50    model->dual();
51#else
52    int numberColumns = model->numberColumns();
53    int numberRows = model->numberRows();
54    // Use dual region
55    double * rhs = model->dualRowSolution();
56    int * whichRow = new int[3*numberRows];
57    int * whichColumn = new int[2*numberColumns];
58    int nBound;
59    ClpSimplex * small = static_cast<ClpSimplexOther *> (model)->crunch(rhs, whichRow, whichColumn,
60                         nBound, false, false);
61    if (small) {
62        small->dual();
63        if (small->problemStatus() == 0) {
64            model->setProblemStatus(0);
65            static_cast<ClpSimplexOther *> (model)->afterCrunch(*small, whichRow, whichColumn, nBound);
66        } else if (small->problemStatus() != 3) {
67            model->setProblemStatus(1);
68        } else {
69            if (small->problemStatus() == 3) {
70                // may be problems
71                small->computeObjectiveValue();
72                model->setObjectiveValue(small->objectiveValue());
73                model->setProblemStatus(3);
74            } else {
75                model->setProblemStatus(3);
76            }
77        }
78        delete small;
79    } else {
80        model->setProblemStatus(1);
81    }
82    delete [] whichRow;
83    delete [] whichColumn;
84#endif
85}
86/*
87  On input
88  doAction - 0 just fix in original and return NULL
89             1 return fixed non-presolved solver
90             2 as one but use presolve Inside this
91             3 use presolve and fix ones with large cost
92             ? do heuristics and set best solution
93             ? do BAB and just set best solution
94             10+ then use lastSolution and relax a few
95             -2 cleanup afterwards if using 2
96  On output - number fixed
97*/
98OsiClpSolverInterface *
99fixVubs(CbcModel & model, int skipZero2,
100        int & doAction,
101        CoinMessageHandler * /*generalMessageHandler*/,
102        const double * lastSolution, double dextra[6],
103        int extra[5])
104{
105    if (doAction == 11 && !lastSolution)
106        lastSolution = model.bestSolution();
107    assert (((doAction >= 0 && doAction <= 3) && !lastSolution) || (doAction == 11 && lastSolution));
108    double fractionIntFixed = dextra[3];
109    double fractionFixed = dextra[4];
110    double fixAbove = dextra[2];
111    double fixAboveValue = (dextra[5] > 0.0) ? dextra[5] : 1.0;
112    double time1 = CoinCpuTime();
113    int leaveIntFree = extra[1];
114    OsiSolverInterface * originalSolver = model.solver();
115    OsiClpSolverInterface * originalClpSolver = dynamic_cast< OsiClpSolverInterface*> (originalSolver);
116    ClpSimplex * originalLpSolver = originalClpSolver->getModelPtr();
117    int * originalColumns = NULL;
118    OsiClpSolverInterface * clpSolver;
119    ClpSimplex * lpSolver;
120    ClpPresolve pinfo;
121    assert(originalSolver->getObjSense() > 0);
122    if (doAction == 2 || doAction == 3) {
123        double * saveLB = NULL;
124        double * saveUB = NULL;
125        int numberColumns = originalLpSolver->numberColumns();
126        if (fixAbove > 0.0) {
127            double time1 = CoinCpuTime();
128            originalClpSolver->initialSolve();
129            printf("first solve took %g seconds\n", CoinCpuTime() - time1);
130            double * columnLower = originalLpSolver->columnLower() ;
131            double * columnUpper = originalLpSolver->columnUpper() ;
132            const double * solution = originalLpSolver->primalColumnSolution();
133            saveLB = CoinCopyOfArray(columnLower, numberColumns);
134            saveUB = CoinCopyOfArray(columnUpper, numberColumns);
135            const double * objective = originalLpSolver->getObjCoefficients() ;
136            int iColumn;
137            int nFix = 0;
138            int nArt = 0;
139            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
140                if (objective[iColumn] > fixAbove) {
141                    if (solution[iColumn] < columnLower[iColumn] + 1.0e-8) {
142                        columnUpper[iColumn] = columnLower[iColumn];
143                        nFix++;
144                    } else {
145                        nArt++;
146                    }
147                } else if (objective[iColumn] < -fixAbove) {
148                    if (solution[iColumn] > columnUpper[iColumn] - 1.0e-8) {
149                        columnLower[iColumn] = columnUpper[iColumn];
150                        nFix++;
151                    } else {
152                        nArt++;
153                    }
154                }
155            }
156            printf("%d artificials fixed, %d left as in solution\n", nFix, nArt);
157            lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e-8, true, 10);
158            if (!lpSolver || doAction == 2) {
159                // take off fixing in original
160                memcpy(columnLower, saveLB, numberColumns*sizeof(double));
161                memcpy(columnUpper, saveUB, numberColumns*sizeof(double));
162            }
163            delete [] saveLB;
164            delete [] saveUB;
165            if (!lpSolver) {
166                // try again
167                pinfo.destroyPresolve();
168                lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e-8, true, 10);
169                assert (lpSolver);
170            }
171        } else {
172            lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e-8, true, 10);
173            assert (lpSolver);
174        }
175        clpSolver = new OsiClpSolverInterface(lpSolver, true);
176        assert(lpSolver == clpSolver->getModelPtr());
177        numberColumns = lpSolver->numberColumns();
178        originalColumns = CoinCopyOfArray(pinfo.originalColumns(), numberColumns);
179        doAction = 1;
180    } else {
181        OsiSolverInterface * solver = originalSolver->clone();
182        clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
183        lpSolver = clpSolver->getModelPtr();
184    }
185    // Tighten bounds
186    lpSolver->tightenPrimalBounds(0.0, 11, true);
187    int numberColumns = clpSolver->getNumCols() ;
188    double * saveColumnLower = CoinCopyOfArray(lpSolver->columnLower(), numberColumns);
189    double * saveColumnUpper = CoinCopyOfArray(lpSolver->columnUpper(), numberColumns);
190    //char generalPrint[200];
191    const double *objective = lpSolver->getObjCoefficients() ;
192    double *columnLower = lpSolver->columnLower() ;
193    double *columnUpper = lpSolver->columnUpper() ;
194    int numberRows = clpSolver->getNumRows();
195    int iRow, iColumn;
196
197    // Row copy
198    CoinPackedMatrix matrixByRow(*clpSolver->getMatrixByRow());
199    const double * elementByRow = matrixByRow.getElements();
200    const int * column = matrixByRow.getIndices();
201    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
202    const int * rowLength = matrixByRow.getVectorLengths();
203
204    // Column copy
205    CoinPackedMatrix  matrixByCol(*clpSolver->getMatrixByCol());
206    //const double * element = matrixByCol.getElements();
207    const int * row = matrixByCol.getIndices();
208    const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
209    const int * columnLength = matrixByCol.getVectorLengths();
210
211    const double * rowLower = clpSolver->getRowLower();
212    const double * rowUpper = clpSolver->getRowUpper();
213
214    // Get maximum size of VUB tree
215    // otherColumn is one fixed to 0 if this one zero
216    int nEl = matrixByCol.getNumElements();
217    CoinBigIndex * fixColumn = new CoinBigIndex [numberColumns+1];
218    int * otherColumn = new int [nEl];
219    int * fix = new int[numberColumns];
220    char * mark = new char [numberColumns];
221    memset(mark, 0, numberColumns);
222    int numberInteger = 0;
223    int numberOther = 0;
224    fixColumn[0] = 0;
225    double large = lpSolver->largeValue(); // treat bounds > this as infinite
226#ifndef NDEBUG
227    double large2 = 1.0e10 * large;
228#endif
229    double tolerance = lpSolver->primalTolerance();
230    int * check = new int[numberRows];
231    for (iRow = 0; iRow < numberRows; iRow++) {
232        check[iRow] = -2; // don't check
233        if (rowLower[iRow] < -1.0e6 && rowUpper[iRow] > 1.0e6)
234            continue;// unlikely
235        // possible row
236        int numberPositive = 0;
237        int iPositive = -1;
238        int numberNegative = 0;
239        int iNegative = -1;
240        CoinBigIndex rStart = rowStart[iRow];
241        CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
242        CoinBigIndex j;
243        int kColumn;
244        for (j = rStart; j < rEnd; ++j) {
245            double value = elementByRow[j];
246            kColumn = column[j];
247            if (columnUpper[kColumn] > columnLower[kColumn]) {
248                if (value > 0.0) {
249                    numberPositive++;
250                    iPositive = kColumn;
251                } else {
252                    numberNegative++;
253                    iNegative = kColumn;
254                }
255            }
256        }
257        if (numberPositive == 1 && numberNegative == 1)
258            check[iRow] = -1; // try both
259        if (numberPositive == 1 && rowLower[iRow] > -1.0e20)
260            check[iRow] = iPositive;
261        else if (numberNegative == 1 && rowUpper[iRow] < 1.0e20)
262            check[iRow] = iNegative;
263    }
264    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
265        fix[iColumn] = -1;
266        if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e-8) {
267            if (clpSolver->isInteger(iColumn))
268                numberInteger++;
269            if (columnLower[iColumn] == 0.0) {
270                bool infeasible = false;
271                fix[iColumn] = 0;
272                // fake upper bound
273                double saveUpper = columnUpper[iColumn];
274                columnUpper[iColumn] = 0.0;
275                for (CoinBigIndex i = columnStart[iColumn];
276                        i < columnStart[iColumn] + columnLength[iColumn]; i++) {
277                    iRow = row[i];
278                    if (check[iRow] != -1 && check[iRow] != iColumn)
279                        continue; // unlikely
280                    // possible row
281                    int infiniteUpper = 0;
282                    int infiniteLower = 0;
283                    double maximumUp = 0.0;
284                    double maximumDown = 0.0;
285                    double newBound;
286                    CoinBigIndex rStart = rowStart[iRow];
287                    CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
288                    CoinBigIndex j;
289                    int kColumn;
290                    // Compute possible lower and upper ranges
291                    for (j = rStart; j < rEnd; ++j) {
292                        double value = elementByRow[j];
293                        kColumn = column[j];
294                        if (value > 0.0) {
295                            if (columnUpper[kColumn] >= large) {
296                                ++infiniteUpper;
297                            } else {
298                                maximumUp += columnUpper[kColumn] * value;
299                            }
300                            if (columnLower[kColumn] <= -large) {
301                                ++infiniteLower;
302                            } else {
303                                maximumDown += columnLower[kColumn] * value;
304                            }
305                        } else if (value < 0.0) {
306                            if (columnUpper[kColumn] >= large) {
307                                ++infiniteLower;
308                            } else {
309                                maximumDown += columnUpper[kColumn] * value;
310                            }
311                            if (columnLower[kColumn] <= -large) {
312                                ++infiniteUpper;
313                            } else {
314                                maximumUp += columnLower[kColumn] * value;
315                            }
316                        }
317                    }
318                    // Build in a margin of error
319                    maximumUp += 1.0e-8 * fabs(maximumUp);
320                    maximumDown -= 1.0e-8 * fabs(maximumDown);
321                    double maxUp = maximumUp + infiniteUpper * 1.0e31;
322                    double maxDown = maximumDown - infiniteLower * 1.0e31;
323                    if (maxUp <= rowUpper[iRow] + tolerance &&
324                            maxDown >= rowLower[iRow] - tolerance) {
325                        //printf("Redundant row in vubs %d\n",iRow);
326                    } else {
327                        if (maxUp < rowLower[iRow] - 100.0*tolerance ||
328                                maxDown > rowUpper[iRow] + 100.0*tolerance) {
329                            infeasible = true;
330                            break;
331                        }
332                        double lower = rowLower[iRow];
333                        double upper = rowUpper[iRow];
334                        for (j = rStart; j < rEnd; ++j) {
335                            double value = elementByRow[j];
336                            kColumn = column[j];
337                            double nowLower = columnLower[kColumn];
338                            double nowUpper = columnUpper[kColumn];
339                            if (value > 0.0) {
340                                // positive value
341                                if (lower > -large) {
342                                    if (!infiniteUpper) {
343                                        assert(nowUpper < large2);
344                                        newBound = nowUpper +
345                                                   (lower - maximumUp) / value;
346                                        // relax if original was large
347                                        if (fabs(maximumUp) > 1.0e8)
348                                            newBound -= 1.0e-12 * fabs(maximumUp);
349                                    } else if (infiniteUpper == 1 && nowUpper > large) {
350                                        newBound = (lower - maximumUp) / value;
351                                        // relax if original was large
352                                        if (fabs(maximumUp) > 1.0e8)
353                                            newBound -= 1.0e-12 * fabs(maximumUp);
354                                    } else {
355                                        newBound = -COIN_DBL_MAX;
356                                    }
357                                    if (newBound > nowLower + 1.0e-12 && newBound > -large) {
358                                        // Tighten the lower bound
359                                        // check infeasible (relaxed)
360                                        if (nowUpper < newBound) {
361                                            if (nowUpper - newBound <
362                                                    -100.0*tolerance) {
363                                                infeasible = true;
364                                                break;
365                                            }
366                                        }
367                                    }
368                                }
369                                if (upper < large) {
370                                    if (!infiniteLower) {
371                                        assert(nowLower > - large2);
372                                        newBound = nowLower +
373                                                   (upper - maximumDown) / value;
374                                        // relax if original was large
375                                        if (fabs(maximumDown) > 1.0e8)
376                                            newBound += 1.0e-12 * fabs(maximumDown);
377                                    } else if (infiniteLower == 1 && nowLower < -large) {
378                                        newBound =   (upper - maximumDown) / value;
379                                        // relax if original was large
380                                        if (fabs(maximumDown) > 1.0e8)
381                                            newBound += 1.0e-12 * fabs(maximumDown);
382                                    } else {
383                                        newBound = COIN_DBL_MAX;
384                                    }
385                                    if (newBound < nowUpper - 1.0e-12 && newBound < large) {
386                                        // Tighten the upper bound
387                                        // check infeasible (relaxed)
388                                        if (nowLower > newBound) {
389                                            if (newBound - nowLower <
390                                                    -100.0*tolerance) {
391                                                infeasible = true;
392                                                break;
393                                            } else {
394                                                newBound = nowLower;
395                                            }
396                                        }
397                                        if (!newBound || (clpSolver->isInteger(kColumn) && newBound < 0.999)) {
398                                            // fix to zero
399                                            if (!mark[kColumn]) {
400                                                otherColumn[numberOther++] = kColumn;
401                                                mark[kColumn] = 1;
402                                                if (check[iRow] == -1)
403                                                    check[iRow] = iColumn;
404                                                else
405                                                    assert(check[iRow] == iColumn);
406                                            }
407                                        }
408                                    }
409                                }
410                            } else {
411                                // negative value
412                                if (lower > -large) {
413                                    if (!infiniteUpper) {
414                                        assert(nowLower < large2);
415                                        newBound = nowLower +
416                                                   (lower - maximumUp) / value;
417                                        // relax if original was large
418                                        if (fabs(maximumUp) > 1.0e8)
419                                            newBound += 1.0e-12 * fabs(maximumUp);
420                                    } else if (infiniteUpper == 1 && nowLower < -large) {
421                                        newBound = (lower - maximumUp) / value;
422                                        // relax if original was large
423                                        if (fabs(maximumUp) > 1.0e8)
424                                            newBound += 1.0e-12 * fabs(maximumUp);
425                                    } else {
426                                        newBound = COIN_DBL_MAX;
427                                    }
428                                    if (newBound < nowUpper - 1.0e-12 && newBound < large) {
429                                        // Tighten the upper bound
430                                        // check infeasible (relaxed)
431                                        if (nowLower > newBound) {
432                                            if (newBound - nowLower <
433                                                    -100.0*tolerance) {
434                                                infeasible = true;
435                                                break;
436                                            } else {
437                                                newBound = nowLower;
438                                            }
439                                        }
440                                        if (!newBound || (clpSolver->isInteger(kColumn) && newBound < 0.999)) {
441                                            // fix to zero
442                                            if (!mark[kColumn]) {
443                                                otherColumn[numberOther++] = kColumn;
444                                                mark[kColumn] = 1;
445                                                if (check[iRow] == -1)
446                                                    check[iRow] = iColumn;
447                                                else
448                                                    assert(check[iRow] == iColumn);
449                                            }
450                                        }
451                                    }
452                                }
453                                if (upper < large) {
454                                    if (!infiniteLower) {
455                                        assert(nowUpper < large2);
456                                        newBound = nowUpper +
457                                                   (upper - maximumDown) / value;
458                                        // relax if original was large
459                                        if (fabs(maximumDown) > 1.0e8)
460                                            newBound -= 1.0e-12 * fabs(maximumDown);
461                                    } else if (infiniteLower == 1 && nowUpper > large) {
462                                        newBound =   (upper - maximumDown) / value;
463                                        // relax if original was large
464                                        if (fabs(maximumDown) > 1.0e8)
465                                            newBound -= 1.0e-12 * fabs(maximumDown);
466                                    } else {
467                                        newBound = -COIN_DBL_MAX;
468                                    }
469                                    if (newBound > nowLower + 1.0e-12 && newBound > -large) {
470                                        // Tighten the lower bound
471                                        // check infeasible (relaxed)
472                                        if (nowUpper < newBound) {
473                                            if (nowUpper - newBound <
474                                                    -100.0*tolerance) {
475                                                infeasible = true;
476                                                break;
477                                            }
478                                        }
479                                    }
480                                }
481                            }
482                        }
483                    }
484                }
485                for (int i = fixColumn[iColumn]; i < numberOther; i++)
486                    mark[otherColumn[i]] = 0;
487                // reset bound unless infeasible
488                if (!infeasible || !clpSolver->isInteger(iColumn))
489                    columnUpper[iColumn] = saveUpper;
490                else if (clpSolver->isInteger(iColumn))
491                    columnLower[iColumn] = 1.0;
492            }
493        }
494        fixColumn[iColumn+1] = numberOther;
495    }
496    delete [] check;
497    delete [] mark;
498    // Now do reverse way
499    int * counts = new int [numberColumns];
500    CoinZeroN(counts, numberColumns);
501    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
502        for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++)
503            counts[otherColumn[i]]++;
504    }
505    numberOther = 0;
506    CoinBigIndex * fixColumn2 = new CoinBigIndex [numberColumns+1];
507    int * otherColumn2 = new int [fixColumn[numberColumns]];
508    fixColumn2[0] = 0;
509    for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
510        numberOther += counts[iColumn];
511        counts[iColumn] = 0;
512        fixColumn2[iColumn+1] = numberOther;
513    }
514    // Create other way
515    for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
516        for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
517            int jColumn = otherColumn[i];
518            CoinBigIndex put = fixColumn2[jColumn] + counts[jColumn];
519            counts[jColumn]++;
520            otherColumn2[put] = iColumn;
521        }
522    }
523    // get top layer i.e. those which are not fixed by any other
524    int kLayer = 0;
525    while (true) {
526        int numberLayered = 0;
527        for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
528            if (fix[iColumn] == kLayer) {
529                for (int i = fixColumn2[iColumn]; i < fixColumn2[iColumn+1]; i++) {
530                    int jColumn = otherColumn2[i];
531                    if (fix[jColumn] == kLayer) {
532                        fix[iColumn] = kLayer + 100;
533                    }
534                }
535            }
536            if (fix[iColumn] == kLayer) {
537                numberLayered++;
538            }
539        }
540        if (numberLayered) {
541            kLayer += 100;
542        } else {
543            break;
544        }
545    }
546    for (int iPass = 0; iPass < 2; iPass++) {
547        for (int jLayer = 0; jLayer < kLayer; jLayer++) {
548            int check[] = { -1, 0, 1, 2, 3, 4, 5, 10, 50, 100, 500, 1000, 5000, 10000, COIN_INT_MAX};
549            int nCheck = static_cast<int> (sizeof(check) / sizeof(int));
550            int countsI[20];
551            int countsC[20];
552            assert (nCheck <= 20);
553            memset(countsI, 0, nCheck*sizeof(int));
554            memset(countsC, 0, nCheck*sizeof(int));
555            check[nCheck-1] = numberColumns;
556            int numberLayered = 0;
557            int numberInteger = 0;
558            for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
559                if (fix[iColumn] == jLayer) {
560                    numberLayered++;
561                    int nFix = fixColumn[iColumn+1] - fixColumn[iColumn];
562                    if (iPass) {
563                        // just integers
564                        nFix = 0;
565                        for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
566                            int jColumn = otherColumn[i];
567                            if (clpSolver->isInteger(jColumn))
568                                nFix++;
569                        }
570                    }
571                    int iFix;
572                    for (iFix = 0; iFix < nCheck; iFix++) {
573                        if (nFix <= check[iFix])
574                            break;
575                    }
576                    assert (iFix < nCheck);
577                    if (clpSolver->isInteger(iColumn)) {
578                        numberInteger++;
579                        countsI[iFix]++;
580                    } else {
581                        countsC[iFix]++;
582                    }
583                }
584            }
585            if (numberLayered) {
586                printf("%d (%d integer) at priority %d\n", numberLayered, numberInteger, 1 + (jLayer / 100));
587                char buffer[50];
588                for (int i = 1; i < nCheck; i++) {
589                    if (countsI[i] || countsC[i]) {
590                        if (i == 1)
591                            sprintf(buffer, " ==    zero            ");
592                        else if (i < nCheck - 1)
593                            sprintf(buffer, "> %6d and <= %6d ", check[i-1], check[i]);
594                        else
595                            sprintf(buffer, "> %6d                ", check[i-1]);
596                        printf("%s %8d integers and %8d continuous\n", buffer, countsI[i], countsC[i]);
597                    }
598                }
599            }
600        }
601    }
602    delete [] counts;
603    // Now do fixing
604    {
605        // switch off presolve and up weight
606        ClpSolve solveOptions;
607        //solveOptions.setPresolveType(ClpSolve::presolveOff,0);
608        solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
609        //solveOptions.setSpecialOption(1,3,30); // sprint
610        int numberColumns = lpSolver->numberColumns();
611        int iColumn;
612        bool allSlack = true;
613        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
614            if (lpSolver->getColumnStatus(iColumn) == ClpSimplex::basic) {
615                allSlack = false;
616                break;
617            }
618        }
619        if (allSlack)
620            solveOptions.setSpecialOption(1, 2, 50); // idiot
621        lpSolver->setInfeasibilityCost(1.0e11);
622        lpSolver->defaultFactorizationFrequency();
623        if (doAction != 11)
624            lpSolver->initialSolve(solveOptions);
625        double * columnLower = lpSolver->columnLower();
626        double * columnUpper = lpSolver->columnUpper();
627        double * fullSolution = lpSolver->primalColumnSolution();
628        const double * dj = lpSolver->dualColumnSolution();
629        int iPass = 0;
630#define MAXPROB 2
631        ClpSimplex models[MAXPROB];
632        int pass[MAXPROB];
633        int kPass = -1;
634        int kLayer = 0;
635        int skipZero = 0;
636        if (skipZero2 == -1)
637            skipZero2 = 40; //-1;
638        /* 0 fixed to 0 by choice
639           1 lb of 1 by choice
640           2 fixed to 0 by another
641           3 as 2 but this go
642           -1 free
643        */
644        char * state = new char [numberColumns];
645        for (iColumn = 0; iColumn < numberColumns; iColumn++)
646            state[iColumn] = -1;
647        while (true) {
648            double largest = -0.1;
649            double smallest = 1.1;
650            int iLargest = -1;
651            int iSmallest = -1;
652            int atZero = 0;
653            int atOne = 0;
654            int toZero = 0;
655            int toOne = 0;
656            int numberFree = 0;
657            int numberGreater = 0;
658            columnLower = lpSolver->columnLower();
659            columnUpper = lpSolver->columnUpper();
660            fullSolution = lpSolver->primalColumnSolution();
661            if (doAction == 11) {
662                {
663                    double * columnLower = lpSolver->columnLower();
664                    double * columnUpper = lpSolver->columnUpper();
665                    //    lpSolver->dual();
666                    memcpy(columnLower, saveColumnLower, numberColumns*sizeof(double));
667                    memcpy(columnUpper, saveColumnUpper, numberColumns*sizeof(double));
668                    //    lpSolver->dual();
669                    int iColumn;
670                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
671                        if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e-8) {
672                            if (clpSolver->isInteger(iColumn)) {
673                                double value = lastSolution[iColumn];
674                                int iValue = static_cast<int> (value + 0.5);
675                                assert (fabs(value - static_cast<double> (iValue)) < 1.0e-3);
676                                assert (iValue >= columnLower[iColumn] &&
677                                        iValue <= columnUpper[iColumn]);
678                                columnLower[iColumn] = iValue;
679                                columnUpper[iColumn] = iValue;
680                            }
681                        }
682                    }
683                    lpSolver->initialSolve(solveOptions);
684                    memcpy(columnLower, saveColumnLower, numberColumns*sizeof(double));
685                    memcpy(columnUpper, saveColumnUpper, numberColumns*sizeof(double));
686                }
687                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
688                    if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e-8) {
689                        if (clpSolver->isInteger(iColumn)) {
690                            double value = lastSolution[iColumn];
691                            int iValue = static_cast<int> (value + 0.5);
692                            assert (fabs(value - static_cast<double> (iValue)) < 1.0e-3);
693                            assert (iValue >= columnLower[iColumn] &&
694                                    iValue <= columnUpper[iColumn]);
695                            if (!fix[iColumn]) {
696                                if (iValue == 0) {
697                                    state[iColumn] = 0;
698                                    assert (!columnLower[iColumn]);
699                                    columnUpper[iColumn] = 0.0;
700                                } else if (iValue == 1) {
701                                    state[iColumn] = 1;
702                                    columnLower[iColumn] = 1.0;
703                                } else {
704                                    // leave fixed
705                                    columnLower[iColumn] = iValue;
706                                    columnUpper[iColumn] = iValue;
707                                }
708                            } else if (iValue == 0) {
709                                state[iColumn] = 10;
710                                columnUpper[iColumn] = 0.0;
711                            } else {
712                                // leave fixed
713                                columnLower[iColumn] = iValue;
714                                columnUpper[iColumn] = iValue;
715                            }
716                        }
717                    }
718                }
719                int jLayer = 0;
720                int nFixed = -1;
721                int nTotalFixed = 0;
722                while (nFixed) {
723                    nFixed = 0;
724                    for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
725                        if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {
726                            for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
727                                int jColumn = otherColumn[i];
728                                if (columnUpper[jColumn]) {
729                                    bool canFix = true;
730                                    for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
731                                        int kColumn = otherColumn2[k];
732                                        if (state[kColumn] == 1) {
733                                            canFix = false;
734                                            break;
735                                        }
736                                    }
737                                    if (canFix) {
738                                        columnUpper[jColumn] = 0.0;
739                                        nFixed++;
740                                    }
741                                }
742                            }
743                        }
744                    }
745                    nTotalFixed += nFixed;
746                    jLayer += 100;
747                }
748                printf("This fixes %d variables in lower priorities\n", nTotalFixed);
749                break;
750            }
751            for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
752                if (!clpSolver->isInteger(iColumn) || fix[iColumn] > kLayer)
753                    continue;
754                // skip if fixes nothing
755                if (fixColumn[iColumn+1] - fixColumn[iColumn] <= skipZero2)
756                    continue;
757                double value = fullSolution[iColumn];
758                if (value > 1.00001) {
759                    numberGreater++;
760                    continue;
761                }
762                double lower = columnLower[iColumn];
763                double upper = columnUpper[iColumn];
764                if (lower == upper) {
765                    if (lower)
766                        atOne++;
767                    else
768                        atZero++;
769                    continue;
770                }
771                if (value < 1.0e-7) {
772                    toZero++;
773                    columnUpper[iColumn] = 0.0;
774                    state[iColumn] = 10;
775                    continue;
776                }
777                if (value > 1.0 - 1.0e-7) {
778                    toOne++;
779                    columnLower[iColumn] = 1.0;
780                    state[iColumn] = 1;
781                    continue;
782                }
783                numberFree++;
784                // skip if fixes nothing
785                if (fixColumn[iColumn+1] - fixColumn[iColumn] <= skipZero)
786                    continue;
787                if (value < smallest) {
788                    smallest = value;
789                    iSmallest = iColumn;
790                }
791                if (value > largest) {
792                    largest = value;
793                    iLargest = iColumn;
794                }
795            }
796            if (toZero || toOne)
797                printf("%d at 0 fixed and %d at one fixed\n", toZero, toOne);
798            printf("%d variables free, %d fixed to 0, %d to 1 - smallest %g, largest %g\n",
799                   numberFree, atZero, atOne, smallest, largest);
800            if (numberGreater && !iPass)
801                printf("%d variables have value > 1.0\n", numberGreater);
802            //skipZero2=0; // leave 0 fixing
803            int jLayer = 0;
804            int nFixed = -1;
805            int nTotalFixed = 0;
806            while (nFixed) {
807                nFixed = 0;
808                for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
809                    if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {
810                        for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
811                            int jColumn = otherColumn[i];
812                            if (columnUpper[jColumn]) {
813                                bool canFix = true;
814                                for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
815                                    int kColumn = otherColumn2[k];
816                                    if (state[kColumn] == 1) {
817                                        canFix = false;
818                                        break;
819                                    }
820                                }
821                                if (canFix) {
822                                    columnUpper[jColumn] = 0.0;
823                                    nFixed++;
824                                }
825                            }
826                        }
827                    }
828                }
829                nTotalFixed += nFixed;
830                jLayer += 100;
831            }
832            printf("This fixes %d variables in lower priorities\n", nTotalFixed);
833            if (iLargest < 0 || numberFree <= leaveIntFree)
834                break;
835            double movement;
836            int way;
837            if (smallest <= 1.0 - largest && smallest < 0.2 && largest < fixAboveValue) {
838                columnUpper[iSmallest] = 0.0;
839                state[iSmallest] = 0;
840                movement = smallest;
841                way = -1;
842            } else {
843                columnLower[iLargest] = 1.0;
844                state[iLargest] = 1;
845                movement = 1.0 - largest;
846                way = 1;
847            }
848            double saveObj = lpSolver->objectiveValue();
849            iPass++;
850            kPass = iPass % MAXPROB;
851            models[kPass] = *lpSolver;
852            if (way == -1) {
853                // fix others
854                for (int i = fixColumn[iSmallest]; i < fixColumn[iSmallest+1]; i++) {
855                    int jColumn = otherColumn[i];
856                    if (state[jColumn] == -1) {
857                        columnUpper[jColumn] = 0.0;
858                        state[jColumn] = 3;
859                    }
860                }
861            }
862            pass[kPass] = iPass;
863            double maxCostUp = COIN_DBL_MAX;
864            objective = lpSolver->getObjCoefficients() ;
865            if (way == -1)
866                maxCostUp = (1.0 - movement) * objective[iSmallest];
867            lpSolver->setDualObjectiveLimit(saveObj + maxCostUp);
868            crunchIt(lpSolver);
869            double moveObj = lpSolver->objectiveValue() - saveObj;
870            printf("movement %s was %g costing %g\n",
871                   (way == -1) ? "down" : "up", movement, moveObj);
872            if (way == -1 && (moveObj >= maxCostUp || lpSolver->status())) {
873                // go up
874                columnLower = models[kPass].columnLower();
875                columnUpper = models[kPass].columnUpper();
876                columnLower[iSmallest] = 1.0;
877                columnUpper[iSmallest] = saveColumnUpper[iSmallest];
878                *lpSolver = models[kPass];
879                columnLower = lpSolver->columnLower();
880                columnUpper = lpSolver->columnUpper();
881                fullSolution = lpSolver->primalColumnSolution();
882                dj = lpSolver->dualColumnSolution();
883                columnLower[iSmallest] = 1.0;
884                columnUpper[iSmallest] = saveColumnUpper[iSmallest];
885                state[iSmallest] = 1;
886                // unfix others
887                for (int i = fixColumn[iSmallest]; i < fixColumn[iSmallest+1]; i++) {
888                    int jColumn = otherColumn[i];
889                    if (state[jColumn] == 3) {
890                        columnUpper[jColumn] = saveColumnUpper[jColumn];
891                        state[jColumn] = -1;
892                    }
893                }
894                crunchIt(lpSolver);
895            }
896            models[kPass] = *lpSolver;
897        }
898        lpSolver->dual();
899        printf("Fixing took %g seconds\n", CoinCpuTime() - time1);
900        columnLower = lpSolver->columnLower();
901        columnUpper = lpSolver->columnUpper();
902        fullSolution = lpSolver->primalColumnSolution();
903        dj = lpSolver->dualColumnSolution();
904        int * sort = new int[numberColumns];
905        double * dsort = new double[numberColumns];
906        int chunk = 20;
907        int iRelax = 0;
908        //double fractionFixed=6.0/8.0;
909        // relax while lots fixed
910        while (true) {
911            if (skipZero2 > 10 && doAction < 10)
912                break;
913            iRelax++;
914            int n = 0;
915            double sum0 = 0.0;
916            double sum00 = 0.0;
917            double sum1 = 0.0;
918            for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
919                if (!clpSolver->isInteger(iColumn) || fix[iColumn] > kLayer)
920                    continue;
921                // skip if fixes nothing
922                if (fixColumn[iColumn+1] - fixColumn[iColumn] == 0 && doAction < 10)
923                    continue;
924                double djValue = dj[iColumn];
925                if (state[iColumn] == 1) {
926                    assert (columnLower[iColumn]);
927                    assert (fullSolution[iColumn] > 0.1);
928                    if (djValue > 0.0) {
929                        //printf("YY dj of %d at %g is %g\n",iColumn,value,djValue);
930                        sum1 += djValue;
931                        sort[n] = iColumn;
932                        dsort[n++] = -djValue;
933                    } else {
934                        //printf("dj of %d at %g is %g\n",iColumn,value,djValue);
935                    }
936                } else if (state[iColumn] == 0 || state[iColumn] == 10) {
937                    assert (fullSolution[iColumn] < 0.1);
938                    assert (!columnUpper[iColumn]);
939                    double otherValue = 0.0;
940                    int nn = 0;
941                    for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
942                        int jColumn = otherColumn[i];
943                        if (columnUpper[jColumn] == 0.0) {
944                            if (dj[jColumn] < -1.0e-5) {
945                                nn++;
946                                otherValue += dj[jColumn]; // really need to look at rest
947                            }
948                        }
949                    }
950                    if (djValue < -1.0e-2 || otherValue < -1.0e-2) {
951                        //printf("XX dj of %d at %g is %g - %d out of %d contribute %g\n",iColumn,value,djValue,
952                        // nn,fixColumn[iColumn+1]-fixColumn[iColumn],otherValue);
953                        if (djValue < 1.0e-8) {
954                            sum0 -= djValue;
955                            sum00 -= otherValue;
956                            sort[n] = iColumn;
957                            if (djValue < -1.0e-2)
958                                dsort[n++] = djValue + otherValue;
959                            else
960                                dsort[n++] = djValue + 0.001 * otherValue;
961                        }
962                    } else {
963                        //printf("dj of %d at %g is %g - no contribution from %d\n",iColumn,value,djValue,
964                        //   fixColumn[iColumn+1]-fixColumn[iColumn]);
965                    }
966                }
967            }
968            CoinSort_2(dsort, dsort + n, sort);
969            double * originalColumnLower = saveColumnLower;
970            double * originalColumnUpper = saveColumnUpper;
971            double * lo = CoinCopyOfArray(columnLower, numberColumns);
972            double * up = CoinCopyOfArray(columnUpper, numberColumns);
973            for (int k = 0; k < CoinMin(chunk, n); k++) {
974                iColumn = sort[k];
975                state[iColumn] = -2;
976            }
977            memcpy(columnLower, originalColumnLower, numberColumns*sizeof(double));
978            memcpy(columnUpper, originalColumnUpper, numberColumns*sizeof(double));
979            int nFixed = 0;
980            int nFixed0 = 0;
981            int nFixed1 = 0;
982            for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
983                if (state[iColumn] == 0 || state[iColumn] == 10) {
984                    columnUpper[iColumn] = 0.0;
985                    assert (lo[iColumn] == 0.0);
986                    nFixed++;
987                    nFixed0++;
988                    for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
989                        int jColumn = otherColumn[i];
990                        if (columnUpper[jColumn]) {
991                            bool canFix = true;
992                            for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
993                                int kColumn = otherColumn2[k];
994                                if (state[kColumn] == 1 || state[kColumn] == -2) {
995                                    canFix = false;
996                                    break;
997                                }
998                            }
999                            if (canFix) {
1000                                columnUpper[jColumn] = 0.0;
1001                                assert (lo[jColumn] == 0.0);
1002                                nFixed++;
1003                            }
1004                        }
1005                    }
1006                } else if (state[iColumn] == 1) {
1007                    columnLower[iColumn] = 1.0;
1008                    nFixed1++;
1009                }
1010            }
1011            printf("%d fixed %d orig 0 %d 1\n", nFixed, nFixed0, nFixed1);
1012            int jLayer = 0;
1013            nFixed = -1;
1014            int nTotalFixed = 0;
1015            while (nFixed) {
1016                nFixed = 0;
1017                for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
1018                    if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {
1019                        for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
1020                            int jColumn = otherColumn[i];
1021                            if (columnUpper[jColumn]) {
1022                                bool canFix = true;
1023                                for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
1024                                    int kColumn = otherColumn2[k];
1025                                    if (state[kColumn] == 1 || state[kColumn] == -2) {
1026                                        canFix = false;
1027                                        break;
1028                                    }
1029                                }
1030                                if (canFix) {
1031                                    columnUpper[jColumn] = 0.0;
1032                                    assert (lo[jColumn] == 0.0);
1033                                    nFixed++;
1034                                }
1035                            }
1036                        }
1037                    }
1038                }
1039                nTotalFixed += nFixed;
1040                jLayer += 100;
1041            }
1042            nFixed = 0;
1043            int nFixedI = 0;
1044            for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
1045                if (columnLower[iColumn] == columnUpper[iColumn]) {
1046                    if (clpSolver->isInteger(iColumn))
1047                        nFixedI++;
1048                    nFixed++;
1049                }
1050            }
1051            printf("This fixes %d variables in lower priorities - total %d (%d integer) - all target %d, int target %d\n",
1052                   nTotalFixed, nFixed, nFixedI, static_cast<int>(fractionFixed*numberColumns), static_cast<int> (fractionIntFixed*numberInteger));
1053            int nBad = 0;
1054            int nRelax = 0;
1055            for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
1056                if (lo[iColumn] < columnLower[iColumn] ||
1057                        up[iColumn] > columnUpper[iColumn]) {
1058                    printf("bad %d old %g %g, new %g %g\n", iColumn, lo[iColumn], up[iColumn],
1059                           columnLower[iColumn], columnUpper[iColumn]);
1060                    nBad++;
1061                }
1062                if (lo[iColumn] > columnLower[iColumn] ||
1063                        up[iColumn] < columnUpper[iColumn]) {
1064                    nRelax++;
1065                }
1066            }
1067            printf("%d relaxed\n", nRelax);
1068            if (iRelax > 20 && nRelax == chunk)
1069                nRelax = 0;
1070            if (iRelax > 50)
1071                nRelax = 0;
1072            assert (!nBad);
1073            delete [] lo;
1074            delete [] up;
1075            lpSolver->primal(1);
1076            if (nFixed < fractionFixed*numberColumns || nFixedI < fractionIntFixed*numberInteger || !nRelax)
1077                break;
1078        }
1079        delete [] state;
1080        delete [] sort;
1081        delete [] dsort;
1082    }
1083    delete [] fix;
1084    delete [] fixColumn;
1085    delete [] otherColumn;
1086    delete [] otherColumn2;
1087    delete [] fixColumn2;
1088    // See if was presolved
1089    if (originalColumns) {
1090        columnLower = lpSolver->columnLower();
1091        columnUpper = lpSolver->columnUpper();
1092        for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
1093            saveColumnLower[iColumn] = columnLower[iColumn];
1094            saveColumnUpper[iColumn] = columnUpper[iColumn];
1095        }
1096        pinfo.postsolve(true);
1097        columnLower = originalLpSolver->columnLower();
1098        columnUpper = originalLpSolver->columnUpper();
1099        double * newColumnLower = lpSolver->columnLower();
1100        double * newColumnUpper = lpSolver->columnUpper();
1101        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1102            int jColumn = originalColumns[iColumn];
1103            columnLower[jColumn] = CoinMax(columnLower[jColumn], newColumnLower[iColumn]);
1104            columnUpper[jColumn] = CoinMin(columnUpper[jColumn], newColumnUpper[iColumn]);
1105        }
1106        numberColumns = originalLpSolver->numberColumns();
1107        delete [] originalColumns;
1108    }
1109    delete [] saveColumnLower;
1110    delete [] saveColumnUpper;
1111    if (!originalColumns) {
1112        // Basis
1113        memcpy(originalLpSolver->statusArray(), lpSolver->statusArray(), numberRows + numberColumns);
1114        memcpy(originalLpSolver->primalColumnSolution(), lpSolver->primalColumnSolution(), numberColumns*sizeof(double));
1115        memcpy(originalLpSolver->primalRowSolution(), lpSolver->primalRowSolution(), numberRows*sizeof(double));
1116        // Fix in solver
1117        columnLower = lpSolver->columnLower();
1118        columnUpper = lpSolver->columnUpper();
1119    }
1120    double * originalColumnLower = originalLpSolver->columnLower();
1121    double * originalColumnUpper = originalLpSolver->columnUpper();
1122    // number fixed
1123    doAction = 0;
1124    for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
1125        originalColumnLower[iColumn] = columnLower[iColumn];
1126        originalColumnUpper[iColumn] = columnUpper[iColumn];
1127        if (columnLower[iColumn] == columnUpper[iColumn])
1128            doAction++;
1129    }
1130    printf("%d fixed by vub preprocessing\n", doAction);
1131    if (originalColumns) {
1132        originalLpSolver->initialSolve();
1133    }
1134    delete clpSolver;
1135    return NULL;
1136}
1137
1138int doHeuristics(CbcModel * model, int type, CbcOrClpParam* parameters_,
1139                 int numberParameters_,int noPrinting_,int initialPumpTune)
1140{
1141#ifdef JJF_ZERO //NEW_STYLE_SOLVER==0
1142    CbcOrClpParam * parameters_ = parameters;
1143    int numberParameters_ = numberParameters;
1144    bool noPrinting_ = noPrinting_;
1145#endif
1146    char generalPrint[10000];
1147    CoinMessages generalMessages = model->messages();
1148    CoinMessageHandler * generalMessageHandler = model->messageHandler();
1149    //generalMessageHandler->setPrefix(false);
1150    bool anyToDo = false;
1151    int logLevel = parameters_[whichParam(CLP_PARAM_INT_LOGLEVEL, numberParameters_, parameters_)].intValue();
1152    int useFpump = parameters_[whichParam(CBC_PARAM_STR_FPUMP, numberParameters_, parameters_)].currentOptionAsInteger();
1153    int useRounding = parameters_[whichParam(CBC_PARAM_STR_ROUNDING, numberParameters_, parameters_)].currentOptionAsInteger();
1154    int useGreedy = parameters_[whichParam(CBC_PARAM_STR_GREEDY, numberParameters_, parameters_)].currentOptionAsInteger();
1155    int useCombine = parameters_[whichParam(CBC_PARAM_STR_COMBINE, numberParameters_, parameters_)].currentOptionAsInteger();
1156    int useCrossover = parameters_[whichParam(CBC_PARAM_STR_CROSSOVER2, numberParameters_, parameters_)].currentOptionAsInteger();
1157    //int usePivotC = parameters_[whichParam(CBC_PARAM_STR_PIVOTANDCOMPLEMENT, numberParameters_, parameters_)].currentOptionAsInteger();
1158    int usePivotF = parameters_[whichParam(CBC_PARAM_STR_PIVOTANDFIX, numberParameters_, parameters_)].currentOptionAsInteger();
1159    int useRand = parameters_[whichParam(CBC_PARAM_STR_RANDROUND, numberParameters_, parameters_)].currentOptionAsInteger();
1160    int useRINS = parameters_[whichParam(CBC_PARAM_STR_RINS, numberParameters_, parameters_)].currentOptionAsInteger();
1161    int useRENS = parameters_[whichParam(CBC_PARAM_STR_RENS, numberParameters_, parameters_)].currentOptionAsInteger();
1162    int useDINS = parameters_[whichParam(CBC_PARAM_STR_DINS, numberParameters_, parameters_)].currentOptionAsInteger();
1163    int useDIVING2 = parameters_[whichParam(CBC_PARAM_STR_DIVINGS, numberParameters_, parameters_)].currentOptionAsInteger();
1164    int useNaive = parameters_[whichParam(CBC_PARAM_STR_NAIVE, numberParameters_, parameters_)].currentOptionAsInteger();
1165    int kType = (type < 10) ? type : 1;
1166    assert (kType == 1 || kType == 2);
1167    // FPump done first as it only works if no solution
1168    if (useFpump >= kType && useFpump <= kType + 1) {
1169        anyToDo = true;
1170        CbcHeuristicFPump heuristic4(*model);
1171        double dextra3 = parameters_[whichParam(CBC_PARAM_DBL_SMALLBAB, numberParameters_, parameters_)].doubleValue();
1172        heuristic4.setFractionSmall(dextra3);
1173        double dextra1 = parameters_[whichParam(CBC_PARAM_DBL_ARTIFICIALCOST, numberParameters_, parameters_)].doubleValue();
1174        if (dextra1)
1175            heuristic4.setArtificialCost(dextra1);
1176        heuristic4.setMaximumPasses(parameters_[whichParam(CBC_PARAM_INT_FPUMPITS, numberParameters_, parameters_)].intValue());
1177        if (parameters_[whichParam(CBC_PARAM_INT_FPUMPITS, numberParameters_, parameters_)].intValue() == 21)
1178            heuristic4.setIterationRatio(1.0);
1179        int pumpTune = parameters_[whichParam(CBC_PARAM_INT_FPUMPTUNE, numberParameters_, parameters_)].intValue();
1180        int pumpTune2 = parameters_[whichParam(CBC_PARAM_INT_FPUMPTUNE2, numberParameters_, parameters_)].intValue();
1181        if (pumpTune > 0) {
1182            bool printStuff = (pumpTune != initialPumpTune || logLevel > 1 || pumpTune2 > 0)
1183                              && !noPrinting_;
1184            if (printStuff) {
1185                generalMessageHandler->message(CBC_GENERAL, generalMessages)
1186                << "Options for feasibility pump - "
1187                << CoinMessageEol;
1188            }
1189            /*
1190            >=10000000 for using obj
1191            >=1000000 use as accumulate switch
1192            >=1000 use index+1 as number of large loops
1193            >=100 use dextra1 as cutoff
1194            %100 == 10,20 etc for experimentation
1195            1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
1196            4 and static continuous, 5 as 3 but no internal integers
1197            6 as 3 but all slack basis!
1198            */
1199            double value = model->solver()->getObjSense() * model->solver()->getObjValue();
1200            int w = pumpTune / 10;
1201            int i = w % 10;
1202            w /= 10;
1203            int c = w % 10;
1204            w /= 10;
1205            int r = w;
1206            int accumulate = r / 1000;
1207            r -= 1000 * accumulate;
1208            if (accumulate >= 10) {
1209                int which = accumulate / 10;
1210                accumulate -= 10 * which;
1211                which--;
1212                // weights and factors
1213                double weight[] = {0.01, 0.01, 0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
1214                double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
1215                heuristic4.setInitialWeight(weight[which]);
1216                heuristic4.setWeightFactor(factor[which]);
1217                if (printStuff) {
1218                    sprintf(generalPrint, "Initial weight for objective %g, decay factor %g",
1219                            weight[which], factor[which]);
1220                    generalMessageHandler->message(CBC_GENERAL, generalMessages)
1221                    << generalPrint
1222                    << CoinMessageEol;
1223                }
1224
1225            }
1226            // fake cutoff
1227            if (c) {
1228                double cutoff;
1229                model->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
1230                cutoff = CoinMin(cutoff, value + 0.05 * fabs(value) * c);
1231                double fakeCutoff = parameters_[whichParam(CBC_PARAM_DBL_FAKECUTOFF, numberParameters_, parameters_)].doubleValue();
1232                if (fakeCutoff)
1233                    cutoff = fakeCutoff;
1234                heuristic4.setFakeCutoff(cutoff);
1235                if (printStuff) {
1236                    sprintf(generalPrint, "Fake cutoff of %g", cutoff);
1237                    generalMessageHandler->message(CBC_GENERAL, generalMessages)
1238                    << generalPrint
1239                    << CoinMessageEol;
1240                }
1241            }
1242            int offRandomEtc = 0;
1243            if (pumpTune2) {
1244                if ((pumpTune2 / 1000) != 0) {
1245                    offRandomEtc = 1000000 * (pumpTune2 / 1000);
1246                    if (printStuff) {
1247                        generalMessageHandler->message(CBC_GENERAL, generalMessages)
1248                        << "Feasibility pump may run twice"
1249                        << CoinMessageEol;
1250                    }
1251                    pumpTune2 = pumpTune2 % 1000;
1252                }
1253                if ((pumpTune2 / 100) != 0) {
1254                    offRandomEtc += 100 * (pumpTune2 / 100);
1255                    if (printStuff) {
1256                        generalMessageHandler->message(CBC_GENERAL, generalMessages)
1257                        << "Not using randomized objective"
1258                        << CoinMessageEol;
1259                    }
1260                }
1261                int maxAllowed = pumpTune2 % 100;
1262                if (maxAllowed) {
1263                    offRandomEtc += 1000 * maxAllowed;
1264                    if (printStuff) {
1265                        sprintf(generalPrint, "Fixing if same for %d passes",
1266                                maxAllowed);
1267                        generalMessageHandler->message(CBC_GENERAL, generalMessages)
1268                        << generalPrint
1269                        << CoinMessageEol;
1270                    }
1271                }
1272            }
1273            if (accumulate) {
1274                heuristic4.setAccumulate(accumulate);
1275                if (printStuff) {
1276                    if (accumulate) {
1277                        sprintf(generalPrint, "Accumulate of %d", accumulate);
1278                        generalMessageHandler->message(CBC_GENERAL, generalMessages)
1279                        << generalPrint
1280                        << CoinMessageEol;
1281                    }
1282                }
1283            }
1284            if (r) {
1285                // also set increment
1286                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
1287                double increment = 0.0;
1288                double fakeIncrement = parameters_[whichParam(CBC_PARAM_DBL_FAKEINCREMENT, numberParameters_, parameters_)].doubleValue();
1289                if (fakeIncrement)
1290                    increment = fakeIncrement;
1291                heuristic4.setAbsoluteIncrement(increment);
1292                heuristic4.setMaximumRetries(r + 1);
1293                if (printStuff) {
1294                    if (increment) {
1295                        sprintf(generalPrint, "Increment of %g", increment);
1296                        generalMessageHandler->message(CBC_GENERAL, generalMessages)
1297                        << generalPrint
1298                        << CoinMessageEol;
1299                    }
1300                    sprintf(generalPrint, "%d retries", r + 1);
1301                    generalMessageHandler->message(CBC_GENERAL, generalMessages)
1302                    << generalPrint
1303                    << CoinMessageEol;
1304                }
1305            }
1306            if (i + offRandomEtc) {
1307                heuristic4.setFeasibilityPumpOptions(i*10 + offRandomEtc);
1308                if (printStuff) {
1309                    sprintf(generalPrint, "Feasibility pump options of %d",
1310                            i*10 + offRandomEtc);
1311                    generalMessageHandler->message(CBC_GENERAL, generalMessages)
1312                    << generalPrint
1313                    << CoinMessageEol;
1314                }
1315            }
1316            pumpTune = pumpTune % 100;
1317            if (pumpTune == 6)
1318                pumpTune = 13;
1319            heuristic4.setWhen((pumpTune % 10) + 10);
1320            if (printStuff) {
1321                sprintf(generalPrint, "Tuning (fixing) %d", pumpTune % 10);
1322                generalMessageHandler->message(CBC_GENERAL, generalMessages)
1323                << generalPrint
1324                << CoinMessageEol;
1325            }
1326        }
1327        heuristic4.setHeuristicName("feasibility pump");
1328        //#define ROLF
1329#ifdef ROLF
1330        CbcHeuristicFPump pump(*model);
1331        pump.setMaximumTime(60);
1332        pump.setMaximumPasses(100);
1333        pump.setMaximumRetries(1);
1334        pump.setFixOnReducedCosts(0);
1335        pump.setHeuristicName("Feasibility pump");
1336        pump.setFractionSmall(1.0);
1337        pump.setWhen(13);
1338        model->addHeuristic(&pump);
1339#else
1340        model->addHeuristic(&heuristic4);
1341#endif
1342    }
1343    if (useRounding >= type && useRounding >= kType && useRounding <= kType + 1) {
1344        CbcRounding heuristic1(*model);
1345        heuristic1.setHeuristicName("rounding");
1346        model->addHeuristic(&heuristic1) ;
1347        anyToDo = true;
1348    }
1349    if (useGreedy >= type && useGreedy >= kType && useGreedy <= kType + 1) {
1350        CbcHeuristicGreedyCover heuristic3(*model);
1351        heuristic3.setHeuristicName("greedy cover");
1352        CbcHeuristicGreedyEquality heuristic3a(*model);
1353        heuristic3a.setHeuristicName("greedy equality");
1354        model->addHeuristic(&heuristic3);
1355        model->addHeuristic(&heuristic3a);
1356        anyToDo = true;
1357    }
1358    if ((useRENS==7 && kType==1) || (useRENS==8 && kType==2)) {
1359        useRENS=1+2*(useRENS-7);
1360        CbcHeuristicRENS heuristic6a(*model);
1361        heuristic6a.setHeuristicName("RENSdj");
1362        heuristic6a.setFractionSmall(0.6/*3.4*/);
1363        heuristic6a.setFeasibilityPumpOptions(3);
1364        heuristic6a.setNumberNodes(10);
1365        heuristic6a.setWhereFrom(4*256+4*1);
1366        heuristic6a.setWhen(2);
1367        heuristic6a.setRensType(1+16);
1368        model->addHeuristic(&heuristic6a) ;
1369        heuristic6a.setHeuristicName("RENSub");
1370        heuristic6a.setFractionSmall(0.4);
1371        heuristic6a.setFeasibilityPumpOptions(1008003);
1372        heuristic6a.setNumberNodes(50);
1373        heuristic6a.setWhereFrom(4*256+4*1);
1374        heuristic6a.setWhen(2);
1375        heuristic6a.setRensType(2+16);
1376        model->addHeuristic(&heuristic6a) ;
1377    }
1378    if (useRENS >= kType && useRENS <= kType + 1) {
1379#ifndef JJF_ONE
1380        CbcHeuristicRENS heuristic6(*model);
1381        heuristic6.setHeuristicName("RENS");
1382        heuristic6.setFractionSmall(0.4);
1383        heuristic6.setFeasibilityPumpOptions(1008003);
1384        int nodes [] = { -2, 50, 50, 50, 200, 1000, 10000};
1385        heuristic6.setNumberNodes(nodes[useRENS]);
1386#else
1387        CbcHeuristicVND heuristic6(*model);
1388        heuristic6.setHeuristicName("VND");
1389        heuristic5.setFractionSmall(0.5);
1390        heuristic5.setDecayFactor(5.0);
1391#endif
1392        model->addHeuristic(&heuristic6) ;
1393        anyToDo = true;
1394    }
1395    if (useNaive >= kType && useNaive <= kType + 1) {
1396        CbcHeuristicNaive heuristic5b(*model);
1397        heuristic5b.setHeuristicName("Naive");
1398        heuristic5b.setFractionSmall(0.4);
1399        heuristic5b.setNumberNodes(50);
1400        model->addHeuristic(&heuristic5b) ;
1401        anyToDo = true;
1402    }
1403    int useDIVING = 0;
1404    {
1405        int useD;
1406        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGV, numberParameters_, parameters_)].currentOptionAsInteger();
1407        useDIVING |= 1 * ((useD >= kType) ? 1 : 0);
1408        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGG, numberParameters_, parameters_)].currentOptionAsInteger();
1409        useDIVING |= 2 * ((useD >= kType) ? 1 : 0);
1410        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGF, numberParameters_, parameters_)].currentOptionAsInteger();
1411        useDIVING |= 4 * ((useD >= kType) ? 1 : 0);
1412        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGC, numberParameters_, parameters_)].currentOptionAsInteger();
1413        useDIVING |= 8 * ((useD >= kType) ? 1 : 0);
1414        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGL, numberParameters_, parameters_)].currentOptionAsInteger();
1415        useDIVING |= 16 * ((useD >= kType) ? 1 : 0);
1416        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGP, numberParameters_, parameters_)].currentOptionAsInteger();
1417        useDIVING |= 32 * ((useD >= kType) ? 1 : 0);
1418    }
1419    if (useDIVING2 >= kType && useDIVING2 <= kType + 1) {
1420        int diveOptions = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();
1421        if (diveOptions < 0 || diveOptions > 10)
1422            diveOptions = 2;
1423        CbcHeuristicJustOne heuristicJustOne(*model);
1424        heuristicJustOne.setHeuristicName("DiveAny");
1425        heuristicJustOne.setWhen(diveOptions);
1426        // add in others
1427        CbcHeuristicDiveCoefficient heuristicDC(*model);
1428        heuristicDC.setHeuristicName("DiveCoefficient");
1429        heuristicJustOne.addHeuristic(&heuristicDC, 1.0) ;
1430        CbcHeuristicDiveFractional heuristicDF(*model);
1431        heuristicDF.setHeuristicName("DiveFractional");
1432        heuristicJustOne.addHeuristic(&heuristicDF, 1.0) ;
1433        CbcHeuristicDiveGuided heuristicDG(*model);
1434        heuristicDG.setHeuristicName("DiveGuided");
1435        heuristicJustOne.addHeuristic(&heuristicDG, 1.0) ;
1436        CbcHeuristicDiveLineSearch heuristicDL(*model);
1437        heuristicDL.setHeuristicName("DiveLineSearch");
1438        heuristicJustOne.addHeuristic(&heuristicDL, 1.0) ;
1439        CbcHeuristicDivePseudoCost heuristicDP(*model);
1440        heuristicDP.setHeuristicName("DivePseudoCost");
1441        heuristicJustOne.addHeuristic(&heuristicDP, 1.0) ;
1442        CbcHeuristicDiveVectorLength heuristicDV(*model);
1443        heuristicDV.setHeuristicName("DiveVectorLength");
1444        heuristicJustOne.addHeuristic(&heuristicDV, 1.0) ;
1445        // Now normalize probabilities
1446        heuristicJustOne.normalizeProbabilities();
1447        model->addHeuristic(&heuristicJustOne) ;
1448    }
1449
1450    if (useDIVING > 0) {
1451        int diveOptions2 = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();
1452        int diveOptions;
1453        if (diveOptions2 > 99) {
1454            // switch on various active set stuff
1455            diveOptions = diveOptions2 % 100;
1456            diveOptions2 -= diveOptions;
1457        } else {
1458            diveOptions = diveOptions2;
1459            diveOptions2 = 0;
1460        }
1461        if (diveOptions < 0 || diveOptions > 9)
1462            diveOptions = 2;
1463        if ((useDIVING&1) != 0) {
1464            CbcHeuristicDiveVectorLength heuristicDV(*model);
1465            heuristicDV.setHeuristicName("DiveVectorLength");
1466            heuristicDV.setWhen(diveOptions);
1467            model->addHeuristic(&heuristicDV) ;
1468        }
1469        if ((useDIVING&2) != 0) {
1470            CbcHeuristicDiveGuided heuristicDG(*model);
1471            heuristicDG.setHeuristicName("DiveGuided");
1472            heuristicDG.setWhen(diveOptions);
1473            model->addHeuristic(&heuristicDG) ;
1474        }
1475        if ((useDIVING&4) != 0) {
1476            CbcHeuristicDiveFractional heuristicDF(*model);
1477            heuristicDF.setHeuristicName("DiveFractional");
1478            heuristicDF.setWhen(diveOptions);
1479            model->addHeuristic(&heuristicDF) ;
1480        }
1481        if ((useDIVING&8) != 0) {
1482            CbcHeuristicDiveCoefficient heuristicDC(*model);
1483            heuristicDC.setHeuristicName("DiveCoefficient");
1484            heuristicDC.setWhen(diveOptions);
1485            model->addHeuristic(&heuristicDC) ;
1486        }
1487        if ((useDIVING&16) != 0) {
1488            CbcHeuristicDiveLineSearch heuristicDL(*model);
1489            heuristicDL.setHeuristicName("DiveLineSearch");
1490            heuristicDL.setWhen(diveOptions);
1491            model->addHeuristic(&heuristicDL) ;
1492        }
1493        if ((useDIVING&32) != 0) {
1494            CbcHeuristicDivePseudoCost heuristicDP(*model);
1495            heuristicDP.setHeuristicName("DivePseudoCost");
1496            heuristicDP.setWhen(diveOptions + diveOptions2);
1497            model->addHeuristic(&heuristicDP) ;
1498        }
1499        anyToDo = true;
1500    }
1501#ifdef JJF_ZERO
1502    if (usePivotC >= type && usePivotC <= kType + 1) {
1503        CbcHeuristicPivotAndComplement heuristic(*model);
1504        heuristic.setHeuristicName("pivot and complement");
1505        heuristic.setFractionSmall(10.0); // normally 0.5
1506        model->addHeuristic(&heuristic);
1507        anyToDo = true;
1508    }
1509#endif
1510    if (usePivotF >= type && usePivotF <= kType + 1) {
1511        CbcHeuristicPivotAndFix heuristic(*model);
1512        heuristic.setHeuristicName("pivot and fix");
1513        heuristic.setFractionSmall(10.0); // normally 0.5
1514        model->addHeuristic(&heuristic);
1515        anyToDo = true;
1516    }
1517    if (useRand >= type && useRand <= kType + 1) {
1518        CbcHeuristicRandRound heuristic(*model);
1519        heuristic.setHeuristicName("randomized rounding");
1520        heuristic.setFractionSmall(10.0); // normally 0.5
1521        model->addHeuristic(&heuristic);
1522        anyToDo = true;
1523    }
1524    if (useDINS >= kType && useDINS <= kType + 1) {
1525        CbcHeuristicDINS heuristic5a(*model);
1526        heuristic5a.setHeuristicName("DINS");
1527        heuristic5a.setFractionSmall(0.6);
1528        if (useDINS < 4)
1529            heuristic5a.setDecayFactor(5.0);
1530        else
1531            heuristic5a.setDecayFactor(1.5);
1532        heuristic5a.setNumberNodes(1000);
1533        model->addHeuristic(&heuristic5a) ;
1534        anyToDo = true;
1535    }
1536    if (useRINS >= kType && useRINS <= kType + 1) {
1537        CbcHeuristicRINS heuristic5(*model);
1538        heuristic5.setHeuristicName("RINS");
1539        if (useRINS < 4) {
1540            heuristic5.setFractionSmall(0.5);
1541            heuristic5.setDecayFactor(5.0);
1542        } else {
1543            heuristic5.setFractionSmall(0.6);
1544            heuristic5.setDecayFactor(1.5);
1545        }
1546        model->addHeuristic(&heuristic5) ;
1547        anyToDo = true;
1548    }
1549    if (useCombine >= kType && useCombine <= kType + 1) {
1550        CbcHeuristicLocal heuristic2(*model);
1551        heuristic2.setHeuristicName("combine solutions");
1552        heuristic2.setFractionSmall(0.5);
1553        heuristic2.setSearchType(1);
1554        model->addHeuristic(&heuristic2);
1555        anyToDo = true;
1556    }
1557    if (useCrossover >= kType && useCrossover <= kType + 1) {
1558        CbcHeuristicCrossover heuristic2a(*model);
1559        heuristic2a.setHeuristicName("crossover");
1560        heuristic2a.setFractionSmall(0.3);
1561        // just fix at lower
1562        heuristic2a.setWhen(11);
1563        model->addHeuristic(&heuristic2a);
1564        model->setMaximumSavedSolutions(5);
1565        anyToDo = true;
1566    }
1567    int heurSwitches = parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() % 100;
1568    if (heurSwitches) {
1569        for (int iHeur = 0; iHeur < model->numberHeuristics(); iHeur++) {
1570            CbcHeuristic * heuristic = model->heuristic(iHeur);
1571            heuristic->setSwitches(heurSwitches);
1572        }
1573    }
1574    if (type == 2 && anyToDo) {
1575        // Do heuristics
1576#ifndef JJF_ONE
1577        // clean copy
1578        CbcModel model2(*model);
1579        // But get rid of heuristics in model
1580        model->doHeuristicsAtRoot(2);
1581        if (logLevel <= 1)
1582            model2.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1583        OsiBabSolver defaultC;
1584        //solver_->setAuxiliaryInfo(&defaultC);
1585        model2.passInSolverCharacteristics(&defaultC);
1586        // Save bounds
1587        int numberColumns = model2.solver()->getNumCols();
1588        model2.createContinuousSolver();
1589        bool cleanModel = !model2.numberIntegers() && !model2.numberObjects();
1590        model2.findIntegers(false);
1591        int heurOptions = (parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() / 100) % 100;
1592        if (heurOptions == 0 || heurOptions == 2) {
1593            model2.doHeuristicsAtRoot(1);
1594        } else if (heurOptions == 1 || heurOptions == 3) {
1595            model2.setMaximumNodes(-1);
1596            CbcStrategyDefault strategy(0, 5, 5);
1597            strategy.setupPreProcessing(1, 0);
1598            model2.setStrategy(strategy);
1599            model2.branchAndBound();
1600        }
1601        if (cleanModel)
1602            model2.zapIntegerInformation(false);
1603        if (model2.bestSolution()) {
1604            double value = model2.getMinimizationObjValue();
1605            model->setCutoff(value);
1606            model->setBestSolution(model2.bestSolution(), numberColumns, value);
1607            model->setSolutionCount(1);
1608            model->setNumberHeuristicSolutions(1);
1609        }
1610#else
1611        if (logLevel <= 1)
1612            model->solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1613        OsiBabSolver defaultC;
1614        //solver_->setAuxiliaryInfo(&defaultC);
1615        model->passInSolverCharacteristics(&defaultC);
1616        // Save bounds
1617        int numberColumns = model->solver()->getNumCols();
1618        model->createContinuousSolver();
1619        bool cleanModel = !model->numberIntegers() && !model->numberObjects();
1620        model->findIntegers(false);
1621        model->doHeuristicsAtRoot(1);
1622        if (cleanModel)
1623            model->zapIntegerInformation(false);
1624#endif
1625        return 0;
1626    } else {
1627        return 0;
1628    }
1629}
1630
Note: See TracBrowser for help on using the repository browser.