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

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

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