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

Last change on this file since 2105 was 2105, checked in by forrest, 4 years ago

mostly reporting options plus a few tweaks

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