source: branches/sandbox/Cbc/src/CbcHeuristicDive.cpp @ 1389

Last change on this file since 1389 was 1315, checked in by forrest, 10 years ago

final changes before cleaning

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#if 0
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#if 1
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#if 1
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}
Note: See TracBrowser for help on using the repository browser.