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

Last change on this file since 1566 was 1432, checked in by bjarni, 10 years ago

Added extra return at end of each source file where needed, to remove possible linefeed conflicts (NightlyBuild? errors)

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