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

Last change on this file since 1499 was 1499, checked in by forrest, 9 years ago

a few more rens heuristics

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