source: stable/2.7/Cbc/src/CbcHeuristicDive.cpp

Last change on this file was 1675, checked in by stefan, 8 years ago

sync with trunk rev1674

File size: 39.1 KB
RevLine 
[1271]1/* $Id: CbcHeuristicDive.cpp 1240 2009-10-02 18:41:44Z forrest $ */
[907]2// Copyright (C) 2008, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1573]4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
[907]6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10
11#include "CbcHeuristicDive.hpp"
12#include "CbcStrategy.hpp"
[1271]13#include "OsiAuxInfo.hpp"
[907]14#include  "CoinTime.hpp"
15
[944]16#ifdef COIN_HAS_CLP
17#include "OsiClpSolverInterface.hpp"
18#endif
19
[916]20//#define DIVE_FIX_BINARY_VARIABLES
[944]21//#define DIVE_DEBUG
[916]22
[907]23// Default Constructor
[1286]24CbcHeuristicDive::CbcHeuristicDive()
25        : CbcHeuristic()
[907]26{
[1286]27    // matrix and row copy will automatically be empty
28    downLocks_ = NULL;
29    upLocks_ = NULL;
30    downArray_ = NULL;
31    upArray_ = NULL;
32    percentageToFix_ = 0.2;
33    maxIterations_ = 100;
34    maxSimplexIterations_ = 10000;
35    maxSimplexIterationsAtRoot_ = 1000000;
36    maxTime_ = 600;
37    whereFrom_ = 255 - 2 - 16 + 256;
[1315]38    decayFactor_ = 1.0;
[907]39}
40
41// Constructor from model
42CbcHeuristicDive::CbcHeuristicDive(CbcModel & model)
[1286]43        : CbcHeuristic(model)
[907]44{
[1286]45    downLocks_ = NULL;
46    upLocks_ = NULL;
47    downArray_ = NULL;
48    upArray_ = NULL;
49    // Get a copy of original matrix
50    assert(model.solver());
51    // model may have empty matrix - wait until setModel
52    const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
53    if (matrix) {
54        matrix_ = *matrix;
55        matrixByRow_ = *model.solver()->getMatrixByRow();
56        validate();
57    }
58    percentageToFix_ = 0.2;
59    maxIterations_ = 100;
60    maxSimplexIterations_ = 10000;
61    maxSimplexIterationsAtRoot_ = 1000000;
62    maxTime_ = 600;
63    whereFrom_ = 255 - 2 - 16 + 256;
[1315]64    decayFactor_ = 1.0;
[907]65}
66
[1286]67// Destructor
[907]68CbcHeuristicDive::~CbcHeuristicDive ()
69{
[1286]70    delete [] downLocks_;
71    delete [] upLocks_;
72    assert (!downArray_);
[907]73}
74
75// Create C++ lines to get to current state
[1286]76void
77CbcHeuristicDive::generateCpp( FILE * fp, const char * heuristic)
[907]78{
[1286]79    // hard coded as CbcHeuristic virtual
80    CbcHeuristic::generateCpp(fp, heuristic);
81    if (percentageToFix_ != 0.2)
82        fprintf(fp, "3  %s.setPercentageToFix(%.f);\n", heuristic, percentageToFix_);
83    else
84        fprintf(fp, "4  %s.setPercentageToFix(%.f);\n", heuristic, percentageToFix_);
85    if (maxIterations_ != 100)
86        fprintf(fp, "3  %s.setMaxIterations(%d);\n", heuristic, maxIterations_);
87    else
88        fprintf(fp, "4  %s.setMaxIterations(%d);\n", heuristic, maxIterations_);
89    if (maxSimplexIterations_ != 10000)
90        fprintf(fp, "3  %s.setMaxSimplexIterations(%d);\n", heuristic, maxSimplexIterations_);
91    else
92        fprintf(fp, "4  %s.setMaxSimplexIterations(%d);\n", heuristic, maxSimplexIterations_);
93    if (maxTime_ != 600)
94        fprintf(fp, "3  %s.setMaxTime(%.2f);\n", heuristic, maxTime_);
95    else
96        fprintf(fp, "4  %s.setMaxTime(%.2f);\n", heuristic, maxTime_);
[907]97}
98
[1286]99// Copy constructor
[907]100CbcHeuristicDive::CbcHeuristicDive(const CbcHeuristicDive & rhs)
[1286]101        :
102        CbcHeuristic(rhs),
103        matrix_(rhs.matrix_),
104        matrixByRow_(rhs.matrixByRow_),
105        percentageToFix_(rhs.percentageToFix_),
106        maxIterations_(rhs.maxIterations_),
107        maxSimplexIterations_(rhs.maxSimplexIterations_),
108        maxSimplexIterationsAtRoot_(rhs.maxSimplexIterationsAtRoot_),
109        maxTime_(rhs.maxTime_)
[907]110{
[1286]111    downArray_ = NULL;
112    upArray_ = NULL;
113    if (rhs.downLocks_) {
114        int numberIntegers = model_->numberIntegers();
115        downLocks_ = CoinCopyOfArray(rhs.downLocks_, numberIntegers);
116        upLocks_ = CoinCopyOfArray(rhs.upLocks_, numberIntegers);
117    } else {
118        downLocks_ = NULL;
119        upLocks_ = NULL;
120    }
[907]121}
122
[1286]123// Assignment operator
124CbcHeuristicDive &
125CbcHeuristicDive::operator=( const CbcHeuristicDive & rhs)
[907]126{
[1286]127    if (this != &rhs) {
128        CbcHeuristic::operator=(rhs);
129        matrix_ = rhs.matrix_;
130        matrixByRow_ = rhs.matrixByRow_;
131        percentageToFix_ = rhs.percentageToFix_;
132        maxIterations_ = rhs.maxIterations_;
133        maxSimplexIterations_ = rhs.maxSimplexIterations_;
134        maxSimplexIterationsAtRoot_ = rhs.maxSimplexIterationsAtRoot_;
135        maxTime_ = rhs.maxTime_;
136        delete [] downLocks_;
137        delete [] upLocks_;
138        if (rhs.downLocks_) {
139            int numberIntegers = model_->numberIntegers();
140            downLocks_ = CoinCopyOfArray(rhs.downLocks_, numberIntegers);
141            upLocks_ = CoinCopyOfArray(rhs.upLocks_, numberIntegers);
142        } else {
143            downLocks_ = NULL;
144            upLocks_ = NULL;
145        }
[907]146    }
[1286]147    return *this;
[907]148}
149
150// Resets stuff if model changes
[1286]151void
[907]152CbcHeuristicDive::resetModel(CbcModel * model)
153{
[1286]154    model_ = model;
155    assert(model_->solver());
156    // Get a copy of original matrix
157    const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
158    // model may have empty matrix - wait until setModel
159    if (matrix) {
160        matrix_ = *matrix;
161        matrixByRow_ = *model->solver()->getMatrixByRow();
162        validate();
163    }
[907]164}
165
166// update model
167void CbcHeuristicDive::setModel(CbcModel * model)
168{
[1286]169    model_ = model;
170    assert(model_->solver());
171    // Get a copy of original matrix
172    const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
173    if (matrix) {
174        matrix_ = *matrix;
175        matrixByRow_ = *model->solver()->getMatrixByRow();
176        // make sure model okay for heuristic
177        validate();
178    }
[907]179}
180
181bool CbcHeuristicDive::canHeuristicRun()
182{
[1286]183    return shouldHeurRun_randomChoice();
[907]184}
185
[916]186inline bool compareBinaryVars(const PseudoReducedCost obj1,
[1286]187                              const PseudoReducedCost obj2)
[916]188{
[1286]189    return obj1.pseudoRedCost > obj2.pseudoRedCost;
[916]190}
191
[961]192// See if diving will give better solution
[907]193// Sets value of solution
194// Returns 1 if solution, 0 if not
195int
196CbcHeuristicDive::solution(double & solutionValue,
[1286]197                           double * betterSolution)
[907]198{
[1675]199    int nodeCount = model_->getNodeCount();
200    if (feasibilityPumpOptions_>0 && (nodeCount % feasibilityPumpOptions_) != 0)
201        return 0;
[944]202#ifdef DIVE_DEBUG
[1286]203    std::cout << "solutionValue = " << solutionValue << std::endl;
[944]204#endif
[1286]205    ++numCouldRun_;
[907]206
[1286]207    // test if the heuristic can run
208    if (!canHeuristicRun())
209        return 0;
[907]210
[1393]211#ifdef JJF_ZERO
[1286]212    // See if to do
213    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
214            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
215        return 0; // switched off
[907]216#endif
217
[944]218#ifdef DIVE_DEBUG
[1286]219    int nRoundInfeasible = 0;
220    int nRoundFeasible = 0;
221    int reasonToStop = 0;
[944]222#endif
[1286]223    double time1 = CoinCpuTime();
224    int numberSimplexIterations = 0;
225    int maxSimplexIterations = (model_->getNodeCount()) ? maxSimplexIterations_
226                               : maxSimplexIterationsAtRoot_;
[907]227
[1286]228    OsiSolverInterface * solver = cloneBut(6); // was model_->solver()->clone();
[1271]229# ifdef COIN_HAS_CLP
[1286]230    OsiClpSolverInterface * clpSolver
[1271]231    = dynamic_cast<OsiClpSolverInterface *> (solver);
[1286]232    if (clpSolver) {
233        // say give up easily
234        ClpSimplex * clpSimplex = clpSolver->getModelPtr();
235        clpSimplex->setMoreSpecialOptions(clpSimplex->moreSpecialOptions() | 64);
236    }
[1271]237# endif
[1286]238    const double * lower = solver->getColLower();
239    const double * upper = solver->getColUpper();
240    const double * rowLower = solver->getRowLower();
241    const double * rowUpper = solver->getRowUpper();
242    const double * solution = solver->getColSolution();
243    const double * objective = solver->getObjCoefficients();
244    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
245    double primalTolerance;
246    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
[907]247
[1286]248    int numberRows = matrix_.getNumRows();
249    assert (numberRows <= solver->getNumRows());
250    int numberIntegers = model_->numberIntegers();
251    const int * integerVariable = model_->integerVariable();
252    double direction = solver->getObjSense(); // 1 for min, -1 for max
253    double newSolutionValue = direction * solver->getObjValue();
254    int returnCode = 0;
255    // Column copy
256    const double * element = matrix_.getElements();
257    const int * row = matrix_.getIndices();
258    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
259    const int * columnLength = matrix_.getVectorLengths();
[1271]260#ifdef DIVE_FIX_BINARY_VARIABLES
[1286]261    // Row copy
262    const double * elementByRow = matrixByRow_.getElements();
263    const int * column = matrixByRow_.getIndices();
264    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
265    const int * rowLength = matrixByRow_.getVectorLengths();
[1271]266#endif
[907]267
[1286]268    // Get solution array for heuristic solution
269    int numberColumns = solver->getNumCols();
270    double * newSolution = new double [numberColumns];
271    memcpy(newSolution, solution, numberColumns*sizeof(double));
[907]272
[1286]273    // vectors to store the latest variables fixed at their bounds
274    int* columnFixed = new int [numberIntegers];
275    double* originalBound = new double [numberIntegers];
276    bool * fixedAtLowerBound = new bool [numberIntegers];
277    PseudoReducedCost * candidate = new PseudoReducedCost [numberIntegers];
278    double * random = new double [numberIntegers];
[907]279
[1286]280    int maxNumberAtBoundToFix = static_cast<int> (floor(percentageToFix_ * numberIntegers));
[907]281
[1286]282    // count how many fractional variables
283    int numberFractionalVariables = 0;
284    for (int i = 0; i < numberIntegers; i++) {
285        random[i] = randomNumberGenerator_.randomDouble() + 0.3;
286        int iColumn = integerVariable[i];
287        double value = newSolution[iColumn];
288        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
289            numberFractionalVariables++;
290        }
[907]291    }
292
[1286]293    const double* reducedCost = NULL;
294    // See if not NLP
295    if (model_->solverCharacteristics()->reducedCostsAccurate())
296        reducedCost = solver->getReducedCost();
[916]297
[1286]298    int iteration = 0;
299    while (numberFractionalVariables) {
300        iteration++;
[907]301
[1286]302        // initialize any data
303        initializeData();
[907]304
[1286]305        // select a fractional variable to bound
306        int bestColumn = -1;
307        int bestRound; // -1 rounds down, +1 rounds up
308        bool canRound = selectVariableToBranch(solver, newSolution,
309                                               bestColumn, bestRound);
310
311        // if the solution is not trivially roundable, we don't try to round;
312        // if the solution is trivially roundable, we try to round. However,
313        // if the rounded solution is worse than the current incumbent,
314        // then we don't round and proceed normally. In this case, the
315        // bestColumn will be a trivially roundable variable
316        if (canRound) {
317            // check if by rounding all fractional variables
318            // we get a solution with an objective value
319            // better than the current best integer solution
320            double delta = 0.0;
321            for (int i = 0; i < numberIntegers; i++) {
322                int iColumn = integerVariable[i];
323                double value = newSolution[iColumn];
324                if (fabs(floor(value + 0.5) - value) > integerTolerance) {
325                    assert(downLocks_[i] == 0 || upLocks_[i] == 0);
326                    double obj = objective[iColumn];
327                    if (downLocks_[i] == 0 && upLocks_[i] == 0) {
328                        if (direction * obj >= 0.0)
329                            delta += (floor(value) - value) * obj;
330                        else
331                            delta += (ceil(value) - value) * obj;
332                    } else if (downLocks_[i] == 0)
333                        delta += (floor(value) - value) * obj;
334                    else
335                        delta += (ceil(value) - value) * obj;
336                }
337            }
338            if (direction*(solver->getObjValue() + delta) < solutionValue) {
[944]339#ifdef DIVE_DEBUG
[1286]340                nRoundFeasible++;
[944]341#endif
[1286]342                // Round all the fractional variables
343                for (int i = 0; i < numberIntegers; i++) {
344                    int iColumn = integerVariable[i];
345                    double value = newSolution[iColumn];
346                    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
347                        assert(downLocks_[i] == 0 || upLocks_[i] == 0);
348                        if (downLocks_[i] == 0 && upLocks_[i] == 0) {
349                            if (direction * objective[iColumn] >= 0.0)
350                                newSolution[iColumn] = floor(value);
351                            else
352                                newSolution[iColumn] = ceil(value);
353                        } else if (downLocks_[i] == 0)
354                            newSolution[iColumn] = floor(value);
355                        else
356                            newSolution[iColumn] = ceil(value);
357                    }
358                }
359                break;
360            }
[944]361#ifdef DIVE_DEBUG
[1286]362            else
363                nRoundInfeasible++;
[944]364#endif
[1286]365        }
[944]366
[1286]367        // do reduced cost fixing
[1103]368#ifdef DIVE_DEBUG
[1286]369        int numberFixed = reducedCostFix(solver);
370        std::cout << "numberReducedCostFixed = " << numberFixed << std::endl;
[1103]371#else
[1286]372        reducedCostFix(solver);
[944]373#endif
[1286]374
375        int numberAtBoundFixed = 0;
[944]376#ifdef DIVE_FIX_BINARY_VARIABLES
[1286]377        // fix binary variables based on pseudo reduced cost
378        if (binVarIndex_.size()) {
379            int cnt = 0;
380            int n = static_cast<int>(binVarIndex_.size());
381            for (int j = 0; j < n; j++) {
382                int iColumn1 = binVarIndex_[j];
383                double value = newSolution[iColumn1];
384                if (fabs(value) <= integerTolerance &&
385                        lower[iColumn1] != upper[iColumn1]) {
386                    double maxPseudoReducedCost = 0.0;
[944]387#ifdef DIVE_DEBUG
[1286]388                    std::cout << "iColumn1 = " << iColumn1 << ", value = " << value << std::endl;
[944]389#endif
[1286]390                    int iRow = vbRowIndex_[j];
391                    double chosenValue = 0.0;
392                    for (int k = rowStart[iRow]; k < rowStart[iRow] + rowLength[iRow]; k++) {
393                        int iColumn2 = column[k];
[944]394#ifdef DIVE_DEBUG
[1286]395                        std::cout << "iColumn2 = " << iColumn2 << std::endl;
[944]396#endif
[1286]397                        if (iColumn1 != iColumn2) {
398                            double pseudoReducedCost = fabs(reducedCost[iColumn2] *
399                                                            elementByRow[k]);
[944]400#ifdef DIVE_DEBUG
[1286]401                            int k2;
402                            for (k2 = rowStart[iRow]; k2 < rowStart[iRow] + rowLength[iRow]; k2++) {
403                                if (column[k2] == iColumn1)
404                                    break;
405                            }
406                            std::cout << "reducedCost[" << iColumn2 << "] = "
407                                      << reducedCost[iColumn2]
408                                      << ", elementByRow[" << iColumn2 << "] = " << elementByRow[k]
409                                      << ", elementByRow[" << iColumn1 << "] = " << elementByRow[k2]
410                                      << ", pseudoRedCost = " << pseudoReducedCost
411                                      << std::endl;
[944]412#endif
[1286]413                            if (pseudoReducedCost > maxPseudoReducedCost)
414                                maxPseudoReducedCost = pseudoReducedCost;
415                        } else {
416                            // save value
417                            chosenValue = fabs(elementByRow[k]);
418                        }
419                    }
420                    assert (chosenValue);
421                    maxPseudoReducedCost /= chosenValue;
[944]422#ifdef DIVE_DEBUG
[1286]423                    std::cout << ", maxPseudoRedCost = " << maxPseudoReducedCost << std::endl;
[944]424#endif
[1286]425                    candidate[cnt].var = iColumn1;
426                    candidate[cnt++].pseudoRedCost = maxPseudoReducedCost;
427                }
428            }
[944]429#ifdef DIVE_DEBUG
[1286]430            std::cout << "candidates for rounding = " << cnt << std::endl;
[944]431#endif
[1286]432            std::sort(candidate, candidate + cnt, compareBinaryVars);
433            for (int i = 0; i < cnt; i++) {
434                int iColumn = candidate[i].var;
435                if (numberAtBoundFixed < maxNumberAtBoundToFix) {
436                    columnFixed[numberAtBoundFixed] = iColumn;
437                    originalBound[numberAtBoundFixed] = upper[iColumn];
438                    fixedAtLowerBound[numberAtBoundFixed] = true;
439                    solver->setColUpper(iColumn, lower[iColumn]);
440                    numberAtBoundFixed++;
441                    if (numberAtBoundFixed == maxNumberAtBoundToFix)
442                        break;
443                }
444            }
445        }
[917]446#endif
[916]447
[1286]448        // fix other integer variables that are at their bounds
449        int cnt = 0;
[1271]450#ifdef GAP
[1286]451        double gap = 1.0e30;
[1271]452#endif
[1286]453        if (reducedCost && true) {
[1393]454#ifndef JJF_ONE
[1286]455            cnt = fixOtherVariables(solver, solution, candidate, random);
[1271]456#else
457#ifdef GAP
[1286]458            double cutoff = model_->getCutoff() ;
459            if (cutoff < 1.0e20 && false) {
460                double direction = solver->getObjSense() ;
461                gap = cutoff - solver->getObjValue() * direction ;
462                gap *= 0.1; // Fix more if plausible
463                double tolerance;
464                solver->getDblParam(OsiDualTolerance, tolerance) ;
465                if (gap <= 0.0)
466                    gap = tolerance;
467                gap += 100.0 * tolerance;
468            }
469            int nOverGap = 0;
[1271]470#endif
[1286]471            int numberFree = 0;
472            int numberFixed = 0;
473            for (int i = 0; i < numberIntegers; i++) {
474                int iColumn = integerVariable[i];
475                if (upper[iColumn] > lower[iColumn]) {
476                    numberFree++;
477                    double value = newSolution[iColumn];
478                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
479                        candidate[cnt].var = iColumn;
480                        candidate[cnt++].pseudoRedCost =
481                            fabs(reducedCost[iColumn] * random[i]);
[1271]482#ifdef GAP
[1286]483                        if (fabs(reducedCost[iColumn]) > gap)
484                            nOverGap++;
[1271]485#endif
[1286]486                    }
487                } else {
488                    numberFixed++;
489                }
490            }
[1271]491#ifdef GAP
[1286]492            int nLeft = maxNumberAtBoundToFix - numberAtBoundFixed;
[1315]493#ifdef CLP_INVESTIGATE4
[1286]494            printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
495                   cutoff, solver->getObjValue(), nOverGap, numberFree, numberFixed);
[1271]496#endif
[1286]497            if (nOverGap > nLeft && true) {
498                nOverGap = CoinMin(nOverGap, nLeft + maxNumberAtBoundToFix / 2);
499                maxNumberAtBoundToFix += nOverGap - nLeft;
500            }
[1271]501#else
[1315]502#ifdef CLP_INVESTIGATE4
[1286]503            printf("cutoff %g obj %g - %d free, %d fixed\n",
504                   model_->getCutoff(), solver->getObjValue(), numberFree, numberFixed);
[1271]505#endif
506#endif
507#endif
[1286]508        } else {
509            for (int i = 0; i < numberIntegers; i++) {
510                int iColumn = integerVariable[i];
511                if (upper[iColumn] > lower[iColumn]) {
512                    double value = newSolution[iColumn];
513                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
514                        candidate[cnt].var = iColumn;
515                        candidate[cnt++].pseudoRedCost = numberIntegers - i;
516                    }
517                }
518            }
519        }
520        std::sort(candidate, candidate + cnt, compareBinaryVars);
521        for (int i = 0; i < cnt; i++) {
522            int iColumn = candidate[i].var;
523            if (upper[iColumn] > lower[iColumn]) {
524                double value = newSolution[iColumn];
525                if (fabs(floor(value + 0.5) - value) <= integerTolerance &&
526                        numberAtBoundFixed < maxNumberAtBoundToFix) {
527                    // fix the variable at one of its bounds
528                    if (fabs(lower[iColumn] - value) <= integerTolerance) {
529                        columnFixed[numberAtBoundFixed] = iColumn;
530                        originalBound[numberAtBoundFixed] = upper[iColumn];
531                        fixedAtLowerBound[numberAtBoundFixed] = true;
532                        solver->setColUpper(iColumn, lower[iColumn]);
533                        numberAtBoundFixed++;
534                    } else if (fabs(upper[iColumn] - value) <= integerTolerance) {
535                        columnFixed[numberAtBoundFixed] = iColumn;
536                        originalBound[numberAtBoundFixed] = lower[iColumn];
537                        fixedAtLowerBound[numberAtBoundFixed] = false;
538                        solver->setColLower(iColumn, upper[iColumn]);
539                        numberAtBoundFixed++;
540                    }
541                    if (numberAtBoundFixed == maxNumberAtBoundToFix)
542                        break;
543                }
544            }
545        }
[944]546#ifdef DIVE_DEBUG
[1286]547        std::cout << "numberAtBoundFixed = " << numberAtBoundFixed << std::endl;
[944]548#endif
[907]549
[1286]550        double originalBoundBestColumn;
551        if (bestColumn >= 0) {
552            if (bestRound < 0) {
553                originalBoundBestColumn = upper[bestColumn];
554                solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
555            } else {
556                originalBoundBestColumn = lower[bestColumn];
557                solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
558            }
559        } else {
560            break;
561        }
562        int originalBestRound = bestRound;
563        int saveModelOptions = model_->specialOptions();
564        while (1) {
[907]565
[1286]566            model_->setSpecialOptions(saveModelOptions | 2048);
567            solver->resolve();
568            model_->setSpecialOptions(saveModelOptions);
569            if (!solver->isAbandoned()) {
570                numberSimplexIterations += solver->getIterationCount();
571            } else {
572                numberSimplexIterations = maxSimplexIterations + 1;
573                break;
574            }
[907]575
[1286]576            if (!solver->isProvenOptimal()) {
577                if (numberAtBoundFixed > 0) {
578                    // Remove the bound fix for variables that were at bounds
579                    for (int i = 0; i < numberAtBoundFixed; i++) {
580                        int iColFixed = columnFixed[i];
581                        if (fixedAtLowerBound[i])
582                            solver->setColUpper(iColFixed, originalBound[i]);
583                        else
584                            solver->setColLower(iColFixed, originalBound[i]);
585                    }
586                    numberAtBoundFixed = 0;
587                } else if (bestRound == originalBestRound) {
588                    bestRound *= (-1);
589                    if (bestRound < 0) {
590                        solver->setColLower(bestColumn, originalBoundBestColumn);
591                        solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
592                    } else {
593                        solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
594                        solver->setColUpper(bestColumn, originalBoundBestColumn);
595                    }
596                } else
597                    break;
598            } else
599                break;
600        }
601
602        if (!solver->isProvenOptimal() ||
603                direction*solver->getObjValue() >= solutionValue) {
[944]604#ifdef DIVE_DEBUG
[1286]605            reasonToStop = 1;
[944]606#endif
[1286]607            break;
608        }
[907]609
[1286]610        if (iteration > maxIterations_) {
[944]611#ifdef DIVE_DEBUG
[1286]612            reasonToStop = 2;
[944]613#endif
[1286]614            break;
615        }
[907]616
[1286]617        if (CoinCpuTime() - time1 > maxTime_) {
[944]618#ifdef DIVE_DEBUG
[1286]619            reasonToStop = 3;
[944]620#endif
[1286]621            break;
622        }
[907]623
[1286]624        if (numberSimplexIterations > maxSimplexIterations) {
[1013]625#ifdef DIVE_DEBUG
[1286]626            reasonToStop = 4;
[1013]627#endif
[1286]628            // also switch off
[1040]629#ifdef CLP_INVESTIGATE
[1286]630            printf("switching off diving as too many iterations %d, %d allowed\n",
631                   numberSimplexIterations, maxSimplexIterations);
[1040]632#endif
[1286]633            when_ = 0;
634            break;
635        }
[1013]636
[1286]637        if (solver->getIterationCount() > 1000 && iteration > 3) {
[1088]638#ifdef DIVE_DEBUG
[1286]639            reasonToStop = 5;
[1088]640#endif
[1286]641            // also switch off
[1088]642#ifdef CLP_INVESTIGATE
[1286]643            printf("switching off diving one iteration took %d iterations (total %d)\n",
644                   solver->getIterationCount(), numberSimplexIterations);
[1088]645#endif
[1286]646            when_ = 0;
647            break;
648        }
[1088]649
[1286]650        memcpy(newSolution, solution, numberColumns*sizeof(double));
651        numberFractionalVariables = 0;
652        for (int i = 0; i < numberIntegers; i++) {
653            int iColumn = integerVariable[i];
654            double value = newSolution[iColumn];
655            if (fabs(floor(value + 0.5) - value) > integerTolerance) {
656                numberFractionalVariables++;
657            }
658        }
659
[907]660    }
661
662
[1286]663    double * rowActivity = new double[numberRows];
664    memset(rowActivity, 0, numberRows*sizeof(double));
[907]665
[1286]666    // re-compute new solution value
667    double objOffset = 0.0;
668    solver->getDblParam(OsiObjOffset, objOffset);
669    newSolutionValue = -objOffset;
670    for (int i = 0 ; i < numberColumns ; i++ )
671        newSolutionValue += objective[i] * newSolution[i];
672    newSolutionValue *= direction;
[907]673    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
[1286]674    if (newSolutionValue < solutionValue) {
675        // paranoid check
676        memset(rowActivity, 0, numberRows*sizeof(double));
677        for (int i = 0; i < numberColumns; i++) {
678            int j;
679            double value = newSolution[i];
680            if (value) {
681                for (j = columnStart[i];
682                        j < columnStart[i] + columnLength[i]; j++) {
683                    int iRow = row[j];
684                    rowActivity[iRow] += value * element[j];
685                }
686            }
687        }
688        // check was approximately feasible
689        bool feasible = true;
690        for (int i = 0; i < numberRows; i++) {
691            if (rowActivity[i] < rowLower[i]) {
692                if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
693                    feasible = false;
694            } else if (rowActivity[i] > rowUpper[i]) {
695                if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
696                    feasible = false;
697            }
698        }
699        for (int i = 0; i < numberIntegers; i++) {
700            int iColumn = integerVariable[i];
701            double value = newSolution[iColumn];
702            if (fabs(floor(value + 0.5) - value) > integerTolerance) {
703                feasible = false;
704                break;
705            }
706        }
707        if (feasible) {
708            // new solution
709            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
710            solutionValue = newSolutionValue;
711            //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
712            returnCode = 1;
713        } else {
714            // Can easily happen
715            //printf("Debug CbcHeuristicDive giving bad solution\n");
716        }
[907]717    }
718
[944]719#ifdef DIVE_DEBUG
[1286]720    std::cout << "nRoundInfeasible = " << nRoundInfeasible
721              << ", nRoundFeasible = " << nRoundFeasible
722              << ", returnCode = " << returnCode
723              << ", reasonToStop = " << reasonToStop
724              << ", simplexIts = " << numberSimplexIterations
725              << ", iterations = " << iteration << std::endl;
[944]726#endif
727
[1286]728    delete [] newSolution;
729    delete [] columnFixed;
730    delete [] originalBound;
731    delete [] fixedAtLowerBound;
732    delete [] candidate;
733    delete [] rowActivity;
734    delete [] random;
735    delete [] downArray_;
736    downArray_ = NULL;
737    delete [] upArray_;
738    upArray_ = NULL;
739    delete solver;
740    return returnCode;
[907]741}
742
743// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
[1286]744void
745CbcHeuristicDive::validate()
[907]746{
[1286]747    if (model_ && (when() % 100) < 10) {
748        if (model_->numberIntegers() !=
749                model_->numberObjects() && (model_->numberObjects() ||
750                                            (model_->specialOptions()&1024) == 0)) {
751            int numberOdd = 0;
752            for (int i = 0; i < model_->numberObjects(); i++) {
753                if (!model_->object(i)->canDoHeuristics())
754                    numberOdd++;
755            }
756            if (numberOdd)
757                setWhen(0);
758        }
[1271]759    }
[907]760
[1286]761    int numberIntegers = model_->numberIntegers();
762    const int * integerVariable = model_->integerVariable();
763    delete [] downLocks_;
764    delete [] upLocks_;
765    downLocks_ = new unsigned short [numberIntegers];
766    upLocks_ = new unsigned short [numberIntegers];
767    // Column copy
768    const double * element = matrix_.getElements();
769    const int * row = matrix_.getIndices();
770    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
771    const int * columnLength = matrix_.getVectorLengths();
772    const double * rowLower = model_->solver()->getRowLower();
773    const double * rowUpper = model_->solver()->getRowUpper();
774    for (int i = 0; i < numberIntegers; i++) {
775        int iColumn = integerVariable[i];
776        int down = 0;
777        int up = 0;
778        if (columnLength[iColumn] > 65535) {
779            setWhen(0);
780            break; // unlikely to work
781        }
782        for (CoinBigIndex j = columnStart[iColumn];
783                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
784            int iRow = row[j];
785            if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
786                up++;
787                down++;
788            } else if (element[j] > 0.0) {
789                if (rowUpper[iRow] < 1.0e20)
790                    up++;
791                else
792                    down++;
793            } else {
794                if (rowLower[iRow] > -1.0e20)
795                    up++;
796                else
797                    down++;
798            }
799        }
800        downLocks_[i] = static_cast<unsigned short> (down);
801        upLocks_[i] = static_cast<unsigned short> (up);
[907]802    }
[916]803
804#ifdef DIVE_FIX_BINARY_VARIABLES
[1286]805    selectBinaryVariables();
[916]806#endif
[907]807}
[916]808
809// Select candidate binary variables for fixing
810void
811CbcHeuristicDive::selectBinaryVariables()
812{
[1286]813    // Row copy
814    const double * elementByRow = matrixByRow_.getElements();
815    const int * column = matrixByRow_.getIndices();
816    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
817    const int * rowLength = matrixByRow_.getVectorLengths();
[916]818
[1286]819    const int numberRows = matrixByRow_.getNumRows();
820    const int numberCols = matrixByRow_.getNumCols();
[916]821
[1286]822    const double * lower = model_->solver()->getColLower();
823    const double * upper = model_->solver()->getColUpper();
824    const double * rowLower = model_->solver()->getRowLower();
825    const double * rowUpper = model_->solver()->getRowUpper();
[916]826
[1286]827    //  const char * integerType = model_->integerType();
[916]828
829
[1286]830    //  const int numberIntegers = model_->numberIntegers();
831    //  const int * integerVariable = model_->integerVariable();
832    const double * objective = model_->solver()->getObjCoefficients();
[916]833
[1286]834    // vector to store the row number of variable bound rows
835    int* rowIndexes = new int [numberCols];
836    memset(rowIndexes, -1, numberCols*sizeof(int));
837
838    for (int i = 0; i < numberRows; i++) {
839        int positiveBinary = -1;
840        int negativeBinary = -1;
841        int nPositiveOther = 0;
842        int nNegativeOther = 0;
843        for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
844            int iColumn = column[k];
845            if (model_->solver()->isInteger(iColumn) &&
846                    lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
847                    objective[iColumn] == 0.0 &&
848                    elementByRow[k] > 0.0 &&
849                    positiveBinary < 0)
850                positiveBinary = iColumn;
851            else if (model_->solver()->isInteger(iColumn) &&
852                     lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
853                     objective[iColumn] == 0.0 &&
854                     elementByRow[k] < 0.0 &&
855                     negativeBinary < 0)
856                negativeBinary = iColumn;
857            else if ((elementByRow[k] > 0.0 &&
858                      lower[iColumn] >= 0.0) ||
859                     (elementByRow[k] < 0.0 &&
860                      upper[iColumn] <= 0.0))
861                nPositiveOther++;
862            else if ((elementByRow[k] > 0.0 &&
863                      lower[iColumn] <= 0.0) ||
864                     (elementByRow[k] < 0.0 &&
865                      upper[iColumn] >= 0.0))
866                nNegativeOther++;
867            if (nPositiveOther > 0 && nNegativeOther > 0)
868                break;
869        }
870        int binVar = -1;
871        if (positiveBinary >= 0 &&
872                (negativeBinary >= 0 || nNegativeOther > 0) &&
873                nPositiveOther == 0 &&
874                rowLower[i] == 0.0 &&
875                rowUpper[i] > 0.0)
876            binVar = positiveBinary;
877        else if (negativeBinary >= 0 &&
878                 (positiveBinary >= 0 || nPositiveOther > 0) &&
879                 nNegativeOther == 0 &&
880                 rowLower[i] < 0.0 &&
881                 rowUpper[i] == 0.0)
882            binVar = negativeBinary;
883        if (binVar >= 0) {
884            if (rowIndexes[binVar] == -1)
885                rowIndexes[binVar] = i;
886            else if (rowIndexes[binVar] >= 0)
887                rowIndexes[binVar] = -2;
888        }
[917]889    }
890
[1286]891    for (int j = 0; j < numberCols; j++) {
892        if (rowIndexes[j] >= 0) {
893            binVarIndex_.push_back(j);
894            vbRowIndex_.push_back(rowIndexes[j]);
895        }
[917]896    }
897
[944]898#ifdef DIVE_DEBUG
[1286]899    std::cout << "number vub Binary = " << binVarIndex_.size() << std::endl;
[944]900#endif
[917]901
[1286]902    delete [] rowIndexes;
903
[917]904}
905
[944]906/*
907  Perform reduced cost fixing on integer variables.
[917]908
[944]909  The variables in question are already nonbasic at bound. We're just nailing
910  down the current situation.
911*/
[917]912
[944]913int CbcHeuristicDive::reducedCostFix (OsiSolverInterface* solver)
[917]914
[944]915{
[1286]916    //return 0; // temp
[1393]917#ifndef JJF_ONE
[1286]918    if (!model_->solverCharacteristics()->reducedCostsAccurate())
919        return 0; //NLP
[944]920#endif
[1286]921    double cutoff = model_->getCutoff() ;
922    if (cutoff > 1.0e20)
923        return 0;
[944]924#ifdef DIVE_DEBUG
[1286]925    std::cout << "cutoff = " << cutoff << std::endl;
[944]926#endif
[1286]927    double direction = solver->getObjSense() ;
928    double gap = cutoff - solver->getObjValue() * direction ;
929    gap *= 0.5; // Fix more
930    double tolerance;
931    solver->getDblParam(OsiDualTolerance, tolerance) ;
932    if (gap <= 0.0)
933        gap = tolerance; //return 0;
934    gap += 100.0 * tolerance;
935    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
[917]936
[1286]937    const double *lower = solver->getColLower() ;
938    const double *upper = solver->getColUpper() ;
939    const double *solution = solver->getColSolution() ;
940    const double *reducedCost = solver->getReducedCost() ;
[917]941
[1286]942    int numberIntegers = model_->numberIntegers();
943    const int * integerVariable = model_->integerVariable();
[917]944
[1286]945    int numberFixed = 0 ;
[917]946
[944]947# ifdef COIN_HAS_CLP
[1286]948    OsiClpSolverInterface * clpSolver
[944]949    = dynamic_cast<OsiClpSolverInterface *> (solver);
[1286]950    ClpSimplex * clpSimplex = NULL;
951    if (clpSolver)
952        clpSimplex = clpSolver->getModelPtr();
[944]953# endif
[1286]954    for (int i = 0 ; i < numberIntegers ; i++) {
955        int iColumn = integerVariable[i] ;
956        double djValue = direction * reducedCost[iColumn] ;
957        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
958            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
[944]959#ifdef COIN_HAS_CLP
[1286]960                // may just have been fixed before
961                if (clpSimplex) {
962                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
[944]963#ifdef COIN_DEVELOP
[1286]964                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
965                               iColumn, clpSimplex->getColumnStatus(iColumn),
966                               djValue, gap, lower[iColumn], upper[iColumn]);
[944]967#endif
[1286]968                    } else {
969                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atLowerBound ||
970                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
971                    }
972                }
[944]973#endif
[1286]974                solver->setColUpper(iColumn, lower[iColumn]) ;
975                numberFixed++ ;
976            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
[944]977#ifdef COIN_HAS_CLP
[1286]978                // may just have been fixed before
979                if (clpSimplex) {
980                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
[944]981#ifdef COIN_DEVELOP
[1286]982                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
983                               iColumn, clpSimplex->getColumnStatus(iColumn),
984                               djValue, gap, lower[iColumn], upper[iColumn]);
[944]985#endif
[1286]986                    } else {
987                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atUpperBound ||
988                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
989                    }
990                }
[944]991#endif
[1286]992                solver->setColLower(iColumn, upper[iColumn]) ;
993                numberFixed++ ;
994            }
995        }
[916]996    }
[1286]997    return numberFixed;
[916]998}
[1271]999// Fix other variables at bounds
[1286]1000int
[1271]1001CbcHeuristicDive::fixOtherVariables(OsiSolverInterface * solver,
[1286]1002                                    const double * solution,
1003                                    PseudoReducedCost * candidate,
1004                                    const double * random)
[1271]1005{
[1286]1006    const double * lower = solver->getColLower();
1007    const double * upper = solver->getColUpper();
1008    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1009    double primalTolerance;
1010    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
[1271]1011
[1286]1012    int numberIntegers = model_->numberIntegers();
1013    const int * integerVariable = model_->integerVariable();
1014    const double* reducedCost = solver->getReducedCost();
1015    // fix other integer variables that are at their bounds
1016    int cnt = 0;
[1271]1017#ifdef GAP
[1286]1018    double direction = solver->getObjSense(); // 1 for min, -1 for max
1019    double gap = 1.0e30;
[1271]1020#endif
1021#ifdef GAP
[1286]1022    double cutoff = model_->getCutoff() ;
1023    if (cutoff < 1.0e20 && false) {
1024        double direction = solver->getObjSense() ;
1025        gap = cutoff - solver->getObjValue() * direction ;
1026        gap *= 0.1; // Fix more if plausible
1027        double tolerance;
1028        solver->getDblParam(OsiDualTolerance, tolerance) ;
1029        if (gap <= 0.0)
1030            gap = tolerance;
1031        gap += 100.0 * tolerance;
1032    }
1033    int nOverGap = 0;
[1271]1034#endif
[1286]1035    int numberFree = 0;
1036    int numberFixedAlready = 0;
1037    for (int i = 0; i < numberIntegers; i++) {
1038        int iColumn = integerVariable[i];
1039        if (upper[iColumn] > lower[iColumn]) {
1040            numberFree++;
1041            double value = solution[iColumn];
1042            if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
1043                candidate[cnt].var = iColumn;
1044                candidate[cnt++].pseudoRedCost =
1045                    fabs(reducedCost[iColumn] * random[i]);
[1271]1046#ifdef GAP
[1286]1047                if (fabs(reducedCost[iColumn]) > gap)
1048                    nOverGap++;
[1271]1049#endif
[1286]1050            }
1051        } else {
1052            numberFixedAlready++;
1053        }
[1271]1054    }
1055#ifdef GAP
[1286]1056    int nLeft = maxNumberToFix - numberFixedAlready;
[1315]1057#ifdef CLP_INVESTIGATE4
[1286]1058    printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
1059           cutoff, solver->getObjValue(), nOverGap, numberFree,
1060           numberFixedAlready);
[1271]1061#endif
[1286]1062    if (nOverGap > nLeft && true) {
1063        nOverGap = CoinMin(nOverGap, nLeft + maxNumberToFix / 2);
1064        maxNumberToFix += nOverGap - nLeft;
1065    }
[1271]1066#else
[1315]1067#ifdef CLP_INVESTIGATE4
[1286]1068    printf("cutoff %g obj %g - %d free, %d fixed\n",
1069           model_->getCutoff(), solver->getObjValue(), numberFree,
1070           numberFixedAlready);
[1271]1071#endif
1072#endif
[1286]1073    return cnt;
[1271]1074}
[1432]1075
Note: See TracBrowser for help on using the repository browser.