source: branches/sandbox/Cbc/src/CbcSolverHeuristics.cpp @ 1386

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

Yet another try at cleaning CbcSolver?. This is an intermediate commit to
reconcile Bjarni's changes with my changes.

File size: 72.6 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#if 0
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#if 0 //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 >= kType && useRENS <= kType + 1) {
1357#if 1
1358        CbcHeuristicRENS heuristic6(*model);
1359        heuristic6.setHeuristicName("RENS");
1360        heuristic6.setFractionSmall(0.4);
1361        heuristic6.setFeasibilityPumpOptions(1008003);
1362        int nodes [] = { -2, 50, 50, 50, 200, 1000, 10000};
1363        heuristic6.setNumberNodes(nodes[useRENS]);
1364#else
1365        CbcHeuristicVND heuristic6(*model);
1366        heuristic6.setHeuristicName("VND");
1367        heuristic5.setFractionSmall(0.5);
1368        heuristic5.setDecayFactor(5.0);
1369#endif
1370        model->addHeuristic(&heuristic6) ;
1371        anyToDo = true;
1372    }
1373    if (useNaive >= kType && useNaive <= kType + 1) {
1374        CbcHeuristicNaive heuristic5b(*model);
1375        heuristic5b.setHeuristicName("Naive");
1376        heuristic5b.setFractionSmall(0.4);
1377        heuristic5b.setNumberNodes(50);
1378        model->addHeuristic(&heuristic5b) ;
1379        anyToDo = true;
1380    }
1381    int useDIVING = 0;
1382    {
1383        int useD;
1384        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGV, numberParameters_, parameters_)].currentOptionAsInteger();
1385        useDIVING |= 1 * ((useD >= kType) ? 1 : 0);
1386        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGG, numberParameters_, parameters_)].currentOptionAsInteger();
1387        useDIVING |= 2 * ((useD >= kType) ? 1 : 0);
1388        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGF, numberParameters_, parameters_)].currentOptionAsInteger();
1389        useDIVING |= 4 * ((useD >= kType) ? 1 : 0);
1390        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGC, numberParameters_, parameters_)].currentOptionAsInteger();
1391        useDIVING |= 8 * ((useD >= kType) ? 1 : 0);
1392        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGL, numberParameters_, parameters_)].currentOptionAsInteger();
1393        useDIVING |= 16 * ((useD >= kType) ? 1 : 0);
1394        useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGP, numberParameters_, parameters_)].currentOptionAsInteger();
1395        useDIVING |= 32 * ((useD >= kType) ? 1 : 0);
1396    }
1397    if (useDIVING2 >= kType && useDIVING2 <= kType + 1) {
1398        int diveOptions = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();
1399        if (diveOptions < 0 || diveOptions > 10)
1400            diveOptions = 2;
1401        CbcHeuristicJustOne heuristicJustOne(*model);
1402        heuristicJustOne.setHeuristicName("DiveAny");
1403        heuristicJustOne.setWhen(diveOptions);
1404        // add in others
1405        CbcHeuristicDiveCoefficient heuristicDC(*model);
1406        heuristicDC.setHeuristicName("DiveCoefficient");
1407        heuristicJustOne.addHeuristic(&heuristicDC, 1.0) ;
1408        CbcHeuristicDiveFractional heuristicDF(*model);
1409        heuristicDF.setHeuristicName("DiveFractional");
1410        heuristicJustOne.addHeuristic(&heuristicDF, 1.0) ;
1411        CbcHeuristicDiveGuided heuristicDG(*model);
1412        heuristicDG.setHeuristicName("DiveGuided");
1413        heuristicJustOne.addHeuristic(&heuristicDG, 1.0) ;
1414        CbcHeuristicDiveLineSearch heuristicDL(*model);
1415        heuristicDL.setHeuristicName("DiveLineSearch");
1416        heuristicJustOne.addHeuristic(&heuristicDL, 1.0) ;
1417        CbcHeuristicDivePseudoCost heuristicDP(*model);
1418        heuristicDP.setHeuristicName("DivePseudoCost");
1419        heuristicJustOne.addHeuristic(&heuristicDP, 1.0) ;
1420        CbcHeuristicDiveVectorLength heuristicDV(*model);
1421        heuristicDV.setHeuristicName("DiveVectorLength");
1422        heuristicJustOne.addHeuristic(&heuristicDV, 1.0) ;
1423        // Now normalize probabilities
1424        heuristicJustOne.normalizeProbabilities();
1425        model->addHeuristic(&heuristicJustOne) ;
1426    }
1427
1428    if (useDIVING > 0) {
1429        int diveOptions2 = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();
1430        int diveOptions;
1431        if (diveOptions2 > 99) {
1432            // switch on various active set stuff
1433            diveOptions = diveOptions2 % 100;
1434            diveOptions2 -= diveOptions;
1435        } else {
1436            diveOptions = diveOptions2;
1437            diveOptions2 = 0;
1438        }
1439        if (diveOptions < 0 || diveOptions > 9)
1440            diveOptions = 2;
1441        if ((useDIVING&1) != 0) {
1442            CbcHeuristicDiveVectorLength heuristicDV(*model);
1443            heuristicDV.setHeuristicName("DiveVectorLength");
1444            heuristicDV.setWhen(diveOptions);
1445            model->addHeuristic(&heuristicDV) ;
1446        }
1447        if ((useDIVING&2) != 0) {
1448            CbcHeuristicDiveGuided heuristicDG(*model);
1449            heuristicDG.setHeuristicName("DiveGuided");
1450            heuristicDG.setWhen(diveOptions);
1451            model->addHeuristic(&heuristicDG) ;
1452        }
1453        if ((useDIVING&4) != 0) {
1454            CbcHeuristicDiveFractional heuristicDF(*model);
1455            heuristicDF.setHeuristicName("DiveFractional");
1456            heuristicDF.setWhen(diveOptions);
1457            model->addHeuristic(&heuristicDF) ;
1458        }
1459        if ((useDIVING&8) != 0) {
1460            CbcHeuristicDiveCoefficient heuristicDC(*model);
1461            heuristicDC.setHeuristicName("DiveCoefficient");
1462            heuristicDC.setWhen(diveOptions);
1463            model->addHeuristic(&heuristicDC) ;
1464        }
1465        if ((useDIVING&16) != 0) {
1466            CbcHeuristicDiveLineSearch heuristicDL(*model);
1467            heuristicDL.setHeuristicName("DiveLineSearch");
1468            heuristicDL.setWhen(diveOptions);
1469            model->addHeuristic(&heuristicDL) ;
1470        }
1471        if ((useDIVING&32) != 0) {
1472            CbcHeuristicDivePseudoCost heuristicDP(*model);
1473            heuristicDP.setHeuristicName("DivePseudoCost");
1474            heuristicDP.setWhen(diveOptions + diveOptions2);
1475            model->addHeuristic(&heuristicDP) ;
1476        }
1477        anyToDo = true;
1478    }
1479#if 0
1480    if (usePivotC >= type && usePivotC <= kType + 1) {
1481        CbcHeuristicPivotAndComplement heuristic(*model);
1482        heuristic.setHeuristicName("pivot and complement");
1483        heuristic.setFractionSmall(10.0); // normally 0.5
1484        model->addHeuristic(&heuristic);
1485        anyToDo = true;
1486    }
1487#endif
1488    if (usePivotF >= type && usePivotF <= kType + 1) {
1489        CbcHeuristicPivotAndFix heuristic(*model);
1490        heuristic.setHeuristicName("pivot and fix");
1491        heuristic.setFractionSmall(10.0); // normally 0.5
1492        model->addHeuristic(&heuristic);
1493        anyToDo = true;
1494    }
1495    if (useRand >= type && useRand <= kType + 1) {
1496        CbcHeuristicRandRound heuristic(*model);
1497        heuristic.setHeuristicName("randomized rounding");
1498        heuristic.setFractionSmall(10.0); // normally 0.5
1499        model->addHeuristic(&heuristic);
1500        anyToDo = true;
1501    }
1502    if (useDINS >= kType && useDINS <= kType + 1) {
1503        CbcHeuristicDINS heuristic5a(*model);
1504        heuristic5a.setHeuristicName("DINS");
1505        heuristic5a.setFractionSmall(0.6);
1506        if (useDINS < 4)
1507            heuristic5a.setDecayFactor(5.0);
1508        else
1509            heuristic5a.setDecayFactor(1.5);
1510        heuristic5a.setNumberNodes(1000);
1511        model->addHeuristic(&heuristic5a) ;
1512        anyToDo = true;
1513    }
1514    if (useRINS >= kType && useRINS <= kType + 1) {
1515        CbcHeuristicRINS heuristic5(*model);
1516        heuristic5.setHeuristicName("RINS");
1517        if (useRINS < 4) {
1518            heuristic5.setFractionSmall(0.5);
1519            heuristic5.setDecayFactor(5.0);
1520        } else {
1521            heuristic5.setFractionSmall(0.6);
1522            heuristic5.setDecayFactor(1.5);
1523        }
1524        model->addHeuristic(&heuristic5) ;
1525        anyToDo = true;
1526    }
1527    if (useCombine >= kType && useCombine <= kType + 1) {
1528        CbcHeuristicLocal heuristic2(*model);
1529        heuristic2.setHeuristicName("combine solutions");
1530        heuristic2.setFractionSmall(0.5);
1531        heuristic2.setSearchType(1);
1532        model->addHeuristic(&heuristic2);
1533        anyToDo = true;
1534    }
1535    if (useCrossover >= kType && useCrossover <= kType + 1) {
1536        CbcHeuristicCrossover heuristic2a(*model);
1537        heuristic2a.setHeuristicName("crossover");
1538        heuristic2a.setFractionSmall(0.3);
1539        // just fix at lower
1540        heuristic2a.setWhen(11);
1541        model->addHeuristic(&heuristic2a);
1542        model->setMaximumSavedSolutions(5);
1543        anyToDo = true;
1544    }
1545    int heurSwitches = parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() % 100;
1546    if (heurSwitches) {
1547        for (int iHeur = 0; iHeur < model->numberHeuristics(); iHeur++) {
1548            CbcHeuristic * heuristic = model->heuristic(iHeur);
1549            heuristic->setSwitches(heurSwitches);
1550        }
1551    }
1552    if (type == 2 && anyToDo) {
1553        // Do heuristics
1554#if 1
1555        // clean copy
1556        CbcModel model2(*model);
1557        // But get rid of heuristics in model
1558        model->doHeuristicsAtRoot(2);
1559        if (logLevel <= 1)
1560            model2.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1561        OsiBabSolver defaultC;
1562        //solver_->setAuxiliaryInfo(&defaultC);
1563        model2.passInSolverCharacteristics(&defaultC);
1564        // Save bounds
1565        int numberColumns = model2.solver()->getNumCols();
1566        model2.createContinuousSolver();
1567        bool cleanModel = !model2.numberIntegers() && !model2.numberObjects();
1568        model2.findIntegers(false);
1569        int heurOptions = (parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() / 100) % 100;
1570        if (heurOptions == 0 || heurOptions == 2) {
1571            model2.doHeuristicsAtRoot(1);
1572        } else if (heurOptions == 1 || heurOptions == 3) {
1573            model2.setMaximumNodes(-1);
1574            CbcStrategyDefault strategy(0, 5, 5);
1575            strategy.setupPreProcessing(1, 0);
1576            model2.setStrategy(strategy);
1577            model2.branchAndBound();
1578        }
1579        if (cleanModel)
1580            model2.zapIntegerInformation(false);
1581        if (model2.bestSolution()) {
1582            double value = model2.getMinimizationObjValue();
1583            model->setCutoff(value);
1584            model->setBestSolution(model2.bestSolution(), numberColumns, value);
1585            model->setSolutionCount(1);
1586            model->setNumberHeuristicSolutions(1);
1587        }
1588#else
1589        if (logLevel <= 1)
1590            model->solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1591        OsiBabSolver defaultC;
1592        //solver_->setAuxiliaryInfo(&defaultC);
1593        model->passInSolverCharacteristics(&defaultC);
1594        // Save bounds
1595        int numberColumns = model->solver()->getNumCols();
1596        model->createContinuousSolver();
1597        bool cleanModel = !model->numberIntegers() && !model->numberObjects();
1598        model->findIntegers(false);
1599        model->doHeuristicsAtRoot(1);
1600        if (cleanModel)
1601            model->zapIntegerInformation(false);
1602#endif
1603        return 0;
1604    } else {
1605        return 0;
1606    }
1607}
1608
Note: See TracBrowser for help on using the repository browser.