source: trunk/Cbc/src/CbcHeuristicDive.cpp @ 1573

Last change on this file since 1573 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

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