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

Last change on this file since 1641 was 1641, checked in by forrest, 8 years ago

out some printf statements

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