source: stable/2.8/Cbc/src/CbcHeuristicDive.cpp @ 1912

Last change on this file since 1912 was 1912, checked in by stefan, 6 years ago

sync with trunk rev 1911

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.6 KB
Line 
1/* $Id: CbcHeuristicDive.cpp 1912 2013-04-26 10:43:35Z stefan $ */
2// Copyright (C) 2008, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10
11#include "CbcHeuristicDive.hpp"
12#include "CbcStrategy.hpp"
13#include "CbcModel.hpp"
14#include "CbcSubProblem.hpp"
15#include "OsiAuxInfo.hpp"
16#include  "CoinTime.hpp"
17
18#ifdef COIN_HAS_CLP
19#include "OsiClpSolverInterface.hpp"
20#endif
21
22//#define DIVE_FIX_BINARY_VARIABLES
23//#define DIVE_DEBUG
24
25// Default Constructor
26CbcHeuristicDive::CbcHeuristicDive()
27        : CbcHeuristic()
28{
29    // matrix and row copy will automatically be empty
30    downLocks_ = NULL;
31    upLocks_ = NULL;
32    downArray_ = NULL;
33    upArray_ = NULL;
34    percentageToFix_ = 0.2;
35    maxIterations_ = 100;
36    maxSimplexIterations_ = 10000;
37    maxSimplexIterationsAtRoot_ = 1000000;
38    maxTime_ = 600;
39    whereFrom_ = 255 - 2 - 16 + 256;
40    decayFactor_ = 1.0;
41}
42
43// Constructor from model
44CbcHeuristicDive::CbcHeuristicDive(CbcModel & model)
45        : CbcHeuristic(model)
46{
47    downLocks_ = NULL;
48    upLocks_ = NULL;
49    downArray_ = NULL;
50    upArray_ = NULL;
51    // Get a copy of original matrix
52    assert(model.solver());
53    // model may have empty matrix - wait until setModel
54    const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
55    if (matrix) {
56        matrix_ = *matrix;
57        matrixByRow_ = *model.solver()->getMatrixByRow();
58        validate();
59    }
60    percentageToFix_ = 0.2;
61    maxIterations_ = 100;
62    maxSimplexIterations_ = 10000;
63    maxSimplexIterationsAtRoot_ = 1000000;
64    maxTime_ = 600;
65    whereFrom_ = 255 - 2 - 16 + 256;
66    decayFactor_ = 1.0;
67}
68
69// Destructor
70CbcHeuristicDive::~CbcHeuristicDive ()
71{
72    delete [] downLocks_;
73    delete [] upLocks_;
74    assert (!downArray_);
75}
76
77// Create C++ lines to get to current state
78void
79CbcHeuristicDive::generateCpp( FILE * fp, const char * heuristic)
80{
81    // hard coded as CbcHeuristic virtual
82    CbcHeuristic::generateCpp(fp, heuristic);
83    if (percentageToFix_ != 0.2)
84        fprintf(fp, "3  %s.setPercentageToFix(%.f);\n", heuristic, percentageToFix_);
85    else
86        fprintf(fp, "4  %s.setPercentageToFix(%.f);\n", heuristic, percentageToFix_);
87    if (maxIterations_ != 100)
88        fprintf(fp, "3  %s.setMaxIterations(%d);\n", heuristic, maxIterations_);
89    else
90        fprintf(fp, "4  %s.setMaxIterations(%d);\n", heuristic, maxIterations_);
91    if (maxSimplexIterations_ != 10000)
92        fprintf(fp, "3  %s.setMaxSimplexIterations(%d);\n", heuristic, maxSimplexIterations_);
93    else
94        fprintf(fp, "4  %s.setMaxSimplexIterations(%d);\n", heuristic, maxSimplexIterations_);
95    if (maxTime_ != 600)
96        fprintf(fp, "3  %s.setMaxTime(%.2f);\n", heuristic, maxTime_);
97    else
98        fprintf(fp, "4  %s.setMaxTime(%.2f);\n", heuristic, maxTime_);
99}
100
101// Copy constructor
102CbcHeuristicDive::CbcHeuristicDive(const CbcHeuristicDive & rhs)
103        :
104        CbcHeuristic(rhs),
105        matrix_(rhs.matrix_),
106        matrixByRow_(rhs.matrixByRow_),
107        percentageToFix_(rhs.percentageToFix_),
108        maxIterations_(rhs.maxIterations_),
109        maxSimplexIterations_(rhs.maxSimplexIterations_),
110        maxSimplexIterationsAtRoot_(rhs.maxSimplexIterationsAtRoot_),
111        maxTime_(rhs.maxTime_)
112{
113    downArray_ = NULL;
114    upArray_ = NULL;
115    if (rhs.downLocks_) {
116        int numberIntegers = model_->numberIntegers();
117        downLocks_ = CoinCopyOfArray(rhs.downLocks_, numberIntegers);
118        upLocks_ = CoinCopyOfArray(rhs.upLocks_, numberIntegers);
119    } else {
120        downLocks_ = NULL;
121        upLocks_ = NULL;
122    }
123}
124
125// Assignment operator
126CbcHeuristicDive &
127CbcHeuristicDive::operator=( const CbcHeuristicDive & rhs)
128{
129    if (this != &rhs) {
130        CbcHeuristic::operator=(rhs);
131        matrix_ = rhs.matrix_;
132        matrixByRow_ = rhs.matrixByRow_;
133        percentageToFix_ = rhs.percentageToFix_;
134        maxIterations_ = rhs.maxIterations_;
135        maxSimplexIterations_ = rhs.maxSimplexIterations_;
136        maxSimplexIterationsAtRoot_ = rhs.maxSimplexIterationsAtRoot_;
137        maxTime_ = rhs.maxTime_;
138        delete [] downLocks_;
139        delete [] upLocks_;
140        if (rhs.downLocks_) {
141            int numberIntegers = model_->numberIntegers();
142            downLocks_ = CoinCopyOfArray(rhs.downLocks_, numberIntegers);
143            upLocks_ = CoinCopyOfArray(rhs.upLocks_, numberIntegers);
144        } else {
145            downLocks_ = NULL;
146            upLocks_ = NULL;
147        }
148    }
149    return *this;
150}
151
152// Resets stuff if model changes
153void
154CbcHeuristicDive::resetModel(CbcModel * model)
155{
156    model_ = model;
157    assert(model_->solver());
158    // Get a copy of original matrix
159    const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
160    // model may have empty matrix - wait until setModel
161    if (matrix) {
162        matrix_ = *matrix;
163        matrixByRow_ = *model->solver()->getMatrixByRow();
164        validate();
165    }
166}
167
168// update model
169void CbcHeuristicDive::setModel(CbcModel * model)
170{
171    model_ = model;
172    assert(model_->solver());
173    // Get a copy of original matrix
174    const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
175    if (matrix) {
176        matrix_ = *matrix;
177        matrixByRow_ = *model->solver()->getMatrixByRow();
178        // make sure model okay for heuristic
179        validate();
180    }
181}
182
183bool CbcHeuristicDive::canHeuristicRun()
184{
185    return shouldHeurRun_randomChoice();
186}
187
188inline bool compareBinaryVars(const PseudoReducedCost obj1,
189                              const PseudoReducedCost obj2)
190{
191    return obj1.pseudoRedCost > obj2.pseudoRedCost;
192}
193
194// inner part of dive
195int 
196CbcHeuristicDive::solution(double & solutionValue, int & numberNodes,
197                           int & numberCuts, OsiRowCut ** cuts,
198                           CbcSubProblem ** & nodes,
199                           double * newSolution)
200{
201#ifdef DIVE_DEBUG
202    int nRoundInfeasible = 0;
203    int nRoundFeasible = 0;
204#endif
205    int reasonToStop = 0;
206    double time1 = CoinCpuTime();
207    int numberSimplexIterations = 0;
208    int maxSimplexIterations = (model_->getNodeCount()) ? maxSimplexIterations_
209                               : maxSimplexIterationsAtRoot_;
210    // but can't be exactly coin_int_max
211    maxSimplexIterations = CoinMin(maxSimplexIterations,COIN_INT_MAX>>3);
212    OsiSolverInterface * solver = cloneBut(6); // was model_->solver()->clone();
213# ifdef COIN_HAS_CLP
214    OsiClpSolverInterface * clpSolver
215    = dynamic_cast<OsiClpSolverInterface *> (solver);
216    if (clpSolver) {
217      ClpSimplex * clpSimplex = clpSolver->getModelPtr();
218      int oneSolveIts = clpSimplex->maximumIterations();
219      oneSolveIts = CoinMin(1000+2*(clpSimplex->numberRows()+clpSimplex->numberColumns()),oneSolveIts);
220      clpSimplex->setMaximumIterations(oneSolveIts);
221      if (!nodes) {
222        // say give up easily
223        clpSimplex->setMoreSpecialOptions(clpSimplex->moreSpecialOptions() | 64);
224      } else {
225        // get ray
226        int specialOptions = clpSimplex->specialOptions();
227        specialOptions &= ~0x3100000;
228        specialOptions |= 32;
229        clpSimplex->setSpecialOptions(specialOptions);
230        clpSolver->setSpecialOptions(clpSolver->specialOptions() | 1048576);
231        if ((model_->moreSpecialOptions()&16777216)!=0) {
232          // cutoff is constraint
233          clpSolver->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
234        }
235      }
236    }
237# endif
238    const double * lower = solver->getColLower();
239    const double * upper = solver->getColUpper();
240    const double * rowLower = solver->getRowLower();
241    const double * rowUpper = solver->getRowUpper();
242    const double * solution = solver->getColSolution();
243    const double * objective = solver->getObjCoefficients();
244    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
245    double primalTolerance;
246    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
247
248    int numberRows = matrix_.getNumRows();
249    assert (numberRows <= solver->getNumRows());
250    int numberIntegers = model_->numberIntegers();
251    const int * integerVariable = model_->integerVariable();
252    double direction = solver->getObjSense(); // 1 for min, -1 for max
253    double newSolutionValue = direction * solver->getObjValue();
254    int returnCode = 0;
255    // Column copy
256    const double * element = matrix_.getElements();
257    const int * row = matrix_.getIndices();
258    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
259    const int * columnLength = matrix_.getVectorLengths();
260#ifdef DIVE_FIX_BINARY_VARIABLES
261    // Row copy
262    const double * elementByRow = matrixByRow_.getElements();
263    const int * column = matrixByRow_.getIndices();
264    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
265    const int * rowLength = matrixByRow_.getVectorLengths();
266#endif
267
268    // Get solution array for heuristic solution
269    int numberColumns = solver->getNumCols();
270    memcpy(newSolution, solution, numberColumns*sizeof(double));
271
272    // vectors to store the latest variables fixed at their bounds
273    int* columnFixed = new int [numberIntegers];
274    double* originalBound = new double [numberIntegers+2*numberColumns];
275    double * lowerBefore = originalBound+numberIntegers;
276    double * upperBefore = lowerBefore+numberColumns;
277    memcpy(lowerBefore,lower,numberColumns*sizeof(double));
278    memcpy(upperBefore,upper,numberColumns*sizeof(double));
279    double * lastDjs=newSolution+numberColumns;
280    bool * fixedAtLowerBound = new bool [numberIntegers];
281    PseudoReducedCost * candidate = new PseudoReducedCost [numberIntegers];
282    double * random = new double [numberIntegers];
283
284    int maxNumberAtBoundToFix = static_cast<int> (floor(percentageToFix_ * numberIntegers));
285    assert (!maxNumberAtBoundToFix||!nodes);
286
287    // count how many fractional variables
288    int numberFractionalVariables = 0;
289    for (int i = 0; i < numberIntegers; i++) {
290        random[i] = randomNumberGenerator_.randomDouble() + 0.3;
291        int iColumn = integerVariable[i];
292        double value = newSolution[iColumn];
293        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
294            numberFractionalVariables++;
295        }
296    }
297
298    const double* reducedCost = NULL;
299    // See if not NLP
300    if (model_->solverCharacteristics()->reducedCostsAccurate())
301        reducedCost = solver->getReducedCost();
302
303    int iteration = 0;
304    while (numberFractionalVariables) {
305        iteration++;
306
307        // initialize any data
308        initializeData();
309
310        // select a fractional variable to bound
311        int bestColumn = -1;
312        int bestRound; // -1 rounds down, +1 rounds up
313        bool canRound = selectVariableToBranch(solver, newSolution,
314                                               bestColumn, bestRound);
315        // if the solution is not trivially roundable, we don't try to round;
316        // if the solution is trivially roundable, we try to round. However,
317        // if the rounded solution is worse than the current incumbent,
318        // then we don't round and proceed normally. In this case, the
319        // bestColumn will be a trivially roundable variable
320        if (canRound) {
321            // check if by rounding all fractional variables
322            // we get a solution with an objective value
323            // better than the current best integer solution
324            double delta = 0.0;
325            for (int i = 0; i < numberIntegers; i++) {
326                int iColumn = integerVariable[i];
327                double value = newSolution[iColumn];
328                if (fabs(floor(value + 0.5) - value) > integerTolerance) {
329                    assert(downLocks_[i] == 0 || upLocks_[i] == 0);
330                    double obj = objective[iColumn];
331                    if (downLocks_[i] == 0 && upLocks_[i] == 0) {
332                        if (direction * obj >= 0.0)
333                            delta += (floor(value) - value) * obj;
334                        else
335                            delta += (ceil(value) - value) * obj;
336                    } else if (downLocks_[i] == 0)
337                        delta += (floor(value) - value) * obj;
338                    else
339                        delta += (ceil(value) - value) * obj;
340                }
341            }
342            if (direction*(solver->getObjValue() + delta) < solutionValue) {
343#ifdef DIVE_DEBUG
344                nRoundFeasible++;
345#endif
346                if (!nodes||bestColumn<0) {
347                  // Round all the fractional variables
348                  for (int i = 0; i < numberIntegers; i++) {
349                    int iColumn = integerVariable[i];
350                    double value = newSolution[iColumn];
351                    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
352                      assert(downLocks_[i] == 0 || upLocks_[i] == 0);
353                      if (downLocks_[i] == 0 && upLocks_[i] == 0) {
354                        if (direction * objective[iColumn] >= 0.0)
355                          newSolution[iColumn] = floor(value);
356                        else
357                          newSolution[iColumn] = ceil(value);
358                      } else if (downLocks_[i] == 0)
359                        newSolution[iColumn] = floor(value);
360                      else
361                        newSolution[iColumn] = ceil(value);
362                    }
363                  }
364                  break;
365                } else {
366                  // can't round if going to use in branching
367                  int i;
368                  for (i = 0; i < numberIntegers; i++) {
369                    int iColumn = integerVariable[i];
370                    double value = newSolution[bestColumn];
371                    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
372                      if (iColumn==bestColumn) {
373                        assert(downLocks_[i] == 0 || upLocks_[i] == 0);
374                        double obj = objective[bestColumn];
375                        if (downLocks_[i] == 0 && upLocks_[i] == 0) {
376                          if (direction * obj >= 0.0)
377                            bestRound=-1;
378                          else
379                            bestRound=1;
380                        } else if (downLocks_[i] == 0)
381                          bestRound=-1;
382                        else
383                          bestRound=1;
384                        break;
385                      }
386                    }
387                  }
388                }
389            }
390#ifdef DIVE_DEBUG
391            else
392                nRoundInfeasible++;
393#endif
394        }
395
396        // do reduced cost fixing
397#ifdef DIVE_DEBUG
398        int numberFixed = reducedCostFix(solver);
399        std::cout << "numberReducedCostFixed = " << numberFixed << std::endl;
400#else
401        reducedCostFix(solver);
402#endif
403
404        int numberAtBoundFixed = 0;
405#ifdef DIVE_FIX_BINARY_VARIABLES
406        // fix binary variables based on pseudo reduced cost
407        if (binVarIndex_.size()) {
408            int cnt = 0;
409            int n = static_cast<int>(binVarIndex_.size());
410            for (int j = 0; j < n; j++) {
411                int iColumn1 = binVarIndex_[j];
412                double value = newSolution[iColumn1];
413                if (fabs(value) <= integerTolerance &&
414                        lower[iColumn1] != upper[iColumn1]) {
415                    double maxPseudoReducedCost = 0.0;
416#ifdef DIVE_DEBUG
417                    std::cout << "iColumn1 = " << iColumn1 << ", value = " << value << std::endl;
418#endif
419                    int iRow = vbRowIndex_[j];
420                    double chosenValue = 0.0;
421                    for (int k = rowStart[iRow]; k < rowStart[iRow] + rowLength[iRow]; k++) {
422                        int iColumn2 = column[k];
423#ifdef DIVE_DEBUG
424                        std::cout << "iColumn2 = " << iColumn2 << std::endl;
425#endif
426                        if (iColumn1 != iColumn2) {
427                            double pseudoReducedCost = fabs(reducedCost[iColumn2] *
428                                                            elementByRow[k]);
429#ifdef DIVE_DEBUG
430                            int k2;
431                            for (k2 = rowStart[iRow]; k2 < rowStart[iRow] + rowLength[iRow]; k2++) {
432                                if (column[k2] == iColumn1)
433                                    break;
434                            }
435                            std::cout << "reducedCost[" << iColumn2 << "] = "
436                                      << reducedCost[iColumn2]
437                                      << ", elementByRow[" << iColumn2 << "] = " << elementByRow[k]
438                                      << ", elementByRow[" << iColumn1 << "] = " << elementByRow[k2]
439                                      << ", pseudoRedCost = " << pseudoReducedCost
440                                      << std::endl;
441#endif
442                            if (pseudoReducedCost > maxPseudoReducedCost)
443                                maxPseudoReducedCost = pseudoReducedCost;
444                        } else {
445                            // save value
446                            chosenValue = fabs(elementByRow[k]);
447                        }
448                    }
449                    assert (chosenValue);
450                    maxPseudoReducedCost /= chosenValue;
451#ifdef DIVE_DEBUG
452                    std::cout << ", maxPseudoRedCost = " << maxPseudoReducedCost << std::endl;
453#endif
454                    candidate[cnt].var = iColumn1;
455                    candidate[cnt++].pseudoRedCost = maxPseudoReducedCost;
456                }
457            }
458#ifdef DIVE_DEBUG
459            std::cout << "candidates for rounding = " << cnt << std::endl;
460#endif
461            std::sort(candidate, candidate + cnt, compareBinaryVars);
462            for (int i = 0; i < cnt; i++) {
463                int iColumn = candidate[i].var;
464                if (numberAtBoundFixed < maxNumberAtBoundToFix) {
465                    columnFixed[numberAtBoundFixed] = iColumn;
466                    originalBound[numberAtBoundFixed] = upper[iColumn];
467                    fixedAtLowerBound[numberAtBoundFixed] = true;
468                    solver->setColUpper(iColumn, lower[iColumn]);
469                    numberAtBoundFixed++;
470                    if (numberAtBoundFixed == maxNumberAtBoundToFix)
471                        break;
472                }
473            }
474        }
475#endif
476
477        // fix other integer variables that are at their bounds
478        int cnt = 0;
479#ifdef GAP
480        double gap = 1.0e30;
481#endif
482        if (reducedCost && true) {
483#ifndef JJF_ONE
484            cnt = fixOtherVariables(solver, solution, candidate, random);
485#else
486#ifdef GAP
487            double cutoff = model_->getCutoff() ;
488            if (cutoff < 1.0e20 && false) {
489                double direction = solver->getObjSense() ;
490                gap = cutoff - solver->getObjValue() * direction ;
491                gap *= 0.1; // Fix more if plausible
492                double tolerance;
493                solver->getDblParam(OsiDualTolerance, tolerance) ;
494                if (gap <= 0.0)
495                    gap = tolerance;
496                gap += 100.0 * tolerance;
497            }
498            int nOverGap = 0;
499#endif
500            int numberFree = 0;
501            int numberFixed = 0;
502            for (int i = 0; i < numberIntegers; i++) {
503                int iColumn = integerVariable[i];
504                if (upper[iColumn] > lower[iColumn]) {
505                    numberFree++;
506                    double value = newSolution[iColumn];
507                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
508                        candidate[cnt].var = iColumn;
509                        candidate[cnt++].pseudoRedCost =
510                            fabs(reducedCost[iColumn] * random[i]);
511#ifdef GAP
512                        if (fabs(reducedCost[iColumn]) > gap)
513                            nOverGap++;
514#endif
515                    }
516                } else {
517                    numberFixed++;
518                }
519            }
520#ifdef GAP
521            int nLeft = maxNumberAtBoundToFix - numberAtBoundFixed;
522#ifdef CLP_INVESTIGATE4
523            printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
524                   cutoff, solver->getObjValue(), nOverGap, numberFree, numberFixed);
525#endif
526            if (nOverGap > nLeft && true) {
527                nOverGap = CoinMin(nOverGap, nLeft + maxNumberAtBoundToFix / 2);
528                maxNumberAtBoundToFix += nOverGap - nLeft;
529            }
530#else
531#ifdef CLP_INVESTIGATE4
532            printf("cutoff %g obj %g - %d free, %d fixed\n",
533                   model_->getCutoff(), solver->getObjValue(), numberFree, numberFixed);
534#endif
535#endif
536#endif
537        } else {
538            for (int i = 0; i < numberIntegers; i++) {
539                int iColumn = integerVariable[i];
540                if (upper[iColumn] > lower[iColumn]) {
541                    double value = newSolution[iColumn];
542                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
543                        candidate[cnt].var = iColumn;
544                        candidate[cnt++].pseudoRedCost = numberIntegers - i;
545                    }
546                }
547            }
548        }
549        std::sort(candidate, candidate + cnt, compareBinaryVars);
550        for (int i = 0; i < cnt; i++) {
551            int iColumn = candidate[i].var;
552            if (upper[iColumn] > lower[iColumn]) {
553                double value = newSolution[iColumn];
554                if (fabs(floor(value + 0.5) - value) <= integerTolerance &&
555                        numberAtBoundFixed < maxNumberAtBoundToFix) {
556                    // fix the variable at one of its bounds
557                    if (fabs(lower[iColumn] - value) <= integerTolerance) {
558                        columnFixed[numberAtBoundFixed] = iColumn;
559                        originalBound[numberAtBoundFixed] = upper[iColumn];
560                        fixedAtLowerBound[numberAtBoundFixed] = true;
561                        solver->setColUpper(iColumn, lower[iColumn]);
562                        numberAtBoundFixed++;
563                    } else if (fabs(upper[iColumn] - value) <= integerTolerance) {
564                        columnFixed[numberAtBoundFixed] = iColumn;
565                        originalBound[numberAtBoundFixed] = lower[iColumn];
566                        fixedAtLowerBound[numberAtBoundFixed] = false;
567                        solver->setColLower(iColumn, upper[iColumn]);
568                        numberAtBoundFixed++;
569                    }
570                    if (numberAtBoundFixed == maxNumberAtBoundToFix)
571                        break;
572                }
573            }
574        }
575#ifdef DIVE_DEBUG
576        std::cout << "numberAtBoundFixed = " << numberAtBoundFixed << std::endl;
577#endif
578
579        double originalBoundBestColumn;
580        double bestColumnValue;
581        int whichWay;
582        if (bestColumn >= 0) {
583            bestColumnValue = newSolution[bestColumn];
584            if (bestRound < 0) {
585                originalBoundBestColumn = upper[bestColumn];
586                solver->setColUpper(bestColumn, floor(bestColumnValue));
587                whichWay=0;
588            } else {
589                originalBoundBestColumn = lower[bestColumn];
590                solver->setColLower(bestColumn, ceil(bestColumnValue));
591                whichWay=1;
592            }
593        } else {
594            break;
595        }
596        int originalBestRound = bestRound;
597        int saveModelOptions = model_->specialOptions();
598       
599        while (1) {
600
601            model_->setSpecialOptions(saveModelOptions | 2048);
602            solver->resolve();
603            model_->setSpecialOptions(saveModelOptions);
604            if (!solver->isAbandoned()&&!solver->isIterationLimitReached()) {
605                numberSimplexIterations += solver->getIterationCount();
606            } else {
607                numberSimplexIterations = maxSimplexIterations + 1;
608                reasonToStop += 100;
609                break;
610            }
611
612            if (!solver->isProvenOptimal()) {
613                if (nodes) {
614                  if (solver->isProvenPrimalInfeasible()) {
615                    if (maxSimplexIterationsAtRoot_!=COIN_INT_MAX) {
616                      // stop now
617                      printf("stopping on first infeasibility\n");
618                      break;
619                    } else if (cuts) {
620                      // can do conflict cut
621                      printf("could do intermediate conflict cut\n");
622                      bool localCut;
623                      OsiRowCut * cut = model_->conflictCut(solver,localCut);
624                      if (cut) {
625                        if (!localCut) {
626                          model_->makePartialCut(cut,solver);
627                          cuts[numberCuts++]=cut;
628                        } else {
629                          delete cut;
630                        }
631                      }
632                    }
633                  } else {
634                    reasonToStop += 10;
635                    break;
636                  }
637                }
638                if (numberAtBoundFixed > 0) {
639                    // Remove the bound fix for variables that were at bounds
640                    for (int i = 0; i < numberAtBoundFixed; i++) {
641                        int iColFixed = columnFixed[i];
642                        if (fixedAtLowerBound[i])
643                            solver->setColUpper(iColFixed, originalBound[i]);
644                        else
645                            solver->setColLower(iColFixed, originalBound[i]);
646                    }
647                    numberAtBoundFixed = 0;
648                } else if (bestRound == originalBestRound) {
649                    bestRound *= (-1);
650                    whichWay |=2;
651                    if (bestRound < 0) {
652                        solver->setColLower(bestColumn, originalBoundBestColumn);
653                        solver->setColUpper(bestColumn, floor(bestColumnValue));
654                    } else {
655                        solver->setColLower(bestColumn, ceil(bestColumnValue));
656                        solver->setColUpper(bestColumn, originalBoundBestColumn);
657                    }
658                } else
659                    break;
660            } else
661                break;
662        }
663
664        if (!solver->isProvenOptimal() ||
665                direction*solver->getObjValue() >= solutionValue) {
666            reasonToStop += 1;
667        } else if (iteration > maxIterations_) {
668            reasonToStop += 2;
669        } else if (CoinCpuTime() - time1 > maxTime_) {
670            reasonToStop += 3;
671        } else if (numberSimplexIterations > maxSimplexIterations) {
672            reasonToStop += 4;
673            // also switch off
674#ifdef CLP_INVESTIGATE
675            printf("switching off diving as too many iterations %d, %d allowed\n",
676                   numberSimplexIterations, maxSimplexIterations);
677#endif
678            when_ = 0;
679        } else if (solver->getIterationCount() > 1000 && iteration > 3 && !nodes) {
680            reasonToStop += 5;
681            // also switch off
682#ifdef CLP_INVESTIGATE
683            printf("switching off diving one iteration took %d iterations (total %d)\n",
684                   solver->getIterationCount(), numberSimplexIterations);
685#endif
686            when_ = 0;
687        }
688
689        memcpy(newSolution, solution, numberColumns*sizeof(double));
690        numberFractionalVariables = 0;
691        double sumFractionalVariables=0.0;
692        for (int i = 0; i < numberIntegers; i++) {
693            int iColumn = integerVariable[i];
694            double value = newSolution[iColumn];
695            double away = fabs(floor(value + 0.5) - value);
696            if (away > integerTolerance) {
697                numberFractionalVariables++;
698                sumFractionalVariables += away;
699            }
700        }
701        if (nodes) {
702          // save information
703          //branchValues[numberNodes]=bestColumnValue;
704          //statuses[numberNodes]=whichWay+(bestColumn<<2);
705          //bases[numberNodes]=solver->getWarmStart();
706          ClpSimplex * simplex = clpSolver->getModelPtr();
707          CbcSubProblem * sub =
708            new CbcSubProblem(clpSolver,lowerBefore,upperBefore,
709                          simplex->statusArray(),numberNodes);
710          nodes[numberNodes]=sub;
711          // other stuff
712          sub->branchValue_=bestColumnValue;
713          sub->problemStatus_=whichWay;
714          sub->branchVariable_=bestColumn;
715          sub->objectiveValue_ = simplex->objectiveValue();
716          sub->sumInfeasibilities_ = sumFractionalVariables;
717          sub->numberInfeasibilities_ = numberFractionalVariables;
718          printf("DiveNode %d column %d way %d bvalue %g obj %g\n",
719                 numberNodes,sub->branchVariable_,sub->problemStatus_,
720                 sub->branchValue_,sub->objectiveValue_);
721          numberNodes++;
722          if (solver->isProvenOptimal()) {
723            memcpy(lastDjs,solver->getReducedCost(),numberColumns*sizeof(double));
724            memcpy(lowerBefore,lower,numberColumns*sizeof(double));
725            memcpy(upperBefore,upper,numberColumns*sizeof(double));
726          }
727        }
728        if (!numberFractionalVariables||reasonToStop)
729          break;
730    }
731    if (nodes) {
732      printf("Exiting dive for reason %d\n",reasonToStop);
733      if (reasonToStop>1) {
734        printf("problems in diving\n");
735        int whichWay=nodes[numberNodes-1]->problemStatus_;
736        CbcSubProblem * sub;
737        if ((whichWay&2)==0) {
738          // leave both ways
739          sub = new CbcSubProblem(*nodes[numberNodes-1]);
740          nodes[numberNodes++]=sub;
741        } else {
742          sub = nodes[numberNodes-1];
743        }
744        if ((whichWay&1)==0)
745          sub->problemStatus_=whichWay|1;
746        else
747          sub->problemStatus_=whichWay&~1;
748      }
749      if (!numberNodes) {
750        // was good at start! - create fake
751        clpSolver->resolve();
752        ClpSimplex * simplex = clpSolver->getModelPtr();
753        CbcSubProblem * sub =
754          new CbcSubProblem(clpSolver,lowerBefore,upperBefore,
755                            simplex->statusArray(),numberNodes);
756        nodes[numberNodes]=sub;
757        // other stuff
758        sub->branchValue_=0.0;
759        sub->problemStatus_=0;
760        sub->branchVariable_=-1;
761        sub->objectiveValue_ = simplex->objectiveValue();
762        sub->sumInfeasibilities_ = 0.0;
763        sub->numberInfeasibilities_ = 0;
764        printf("DiveNode %d column %d way %d bvalue %g obj %g\n",
765               numberNodes,sub->branchVariable_,sub->problemStatus_,
766               sub->branchValue_,sub->objectiveValue_);
767        numberNodes++;
768        assert (solver->isProvenOptimal());
769      }
770      nodes[numberNodes-1]->problemStatus_ |= 256*reasonToStop;
771      // use djs as well
772      if (solver->isProvenPrimalInfeasible()&&cuts) {
773        // can do conflict cut and re-order
774        printf("could do final conflict cut\n");
775        bool localCut;
776        OsiRowCut * cut = model_->conflictCut(solver,localCut);
777        if (cut) {
778          printf("cut - need to use conflict and previous djs\n");
779          if (!localCut) {
780            model_->makePartialCut(cut,solver);
781            cuts[numberCuts++]=cut;
782          } else {
783            delete cut;
784          }
785        } else {
786          printf("bad conflict - just use previous djs\n");
787        }
788      }
789    }
790   
791    // re-compute new solution value
792    double objOffset = 0.0;
793    solver->getDblParam(OsiObjOffset, objOffset);
794    newSolutionValue = -objOffset;
795    for (int i = 0 ; i < numberColumns ; i++ )
796      newSolutionValue += objective[i] * newSolution[i];
797    newSolutionValue *= direction;
798    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
799    if (newSolutionValue < solutionValue && !reasonToStop) {
800      double * rowActivity = new double[numberRows];
801      memset(rowActivity, 0, numberRows*sizeof(double));
802      // paranoid check
803      memset(rowActivity, 0, numberRows*sizeof(double));
804      for (int i = 0; i < numberColumns; i++) {
805        int j;
806        double value = newSolution[i];
807        if (value) {
808          for (j = columnStart[i];
809               j < columnStart[i] + columnLength[i]; j++) {
810            int iRow = row[j];
811            rowActivity[iRow] += value * element[j];
812          }
813        }
814      }
815      // check was approximately feasible
816      bool feasible = true;
817      for (int i = 0; i < numberRows; i++) {
818        if (rowActivity[i] < rowLower[i]) {
819          if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
820            feasible = false;
821        } else if (rowActivity[i] > rowUpper[i]) {
822          if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
823            feasible = false;
824        }
825      }
826      for (int i = 0; i < numberIntegers; i++) {
827        int iColumn = integerVariable[i];
828        double value = newSolution[iColumn];
829        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
830          feasible = false;
831          break;
832        }
833      }
834      if (feasible) {
835        // new solution
836        solutionValue = newSolutionValue;
837        //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
838        //if (cuts)
839        //clpSolver->getModelPtr()->writeMps("good8.mps", 2);
840        returnCode = 1;
841      } else {
842        // Can easily happen
843        //printf("Debug CbcHeuristicDive giving bad solution\n");
844      }
845      delete [] rowActivity;
846    }
847
848#ifdef DIVE_DEBUG
849    std::cout << "nRoundInfeasible = " << nRoundInfeasible
850              << ", nRoundFeasible = " << nRoundFeasible
851              << ", returnCode = " << returnCode
852              << ", reasonToStop = " << reasonToStop
853              << ", simplexIts = " << numberSimplexIterations
854              << ", iterations = " << iteration << std::endl;
855#endif
856
857    delete [] columnFixed;
858    delete [] originalBound;
859    delete [] fixedAtLowerBound;
860    delete [] candidate;
861    delete [] random;
862    delete [] downArray_;
863    downArray_ = NULL;
864    delete [] upArray_;
865    upArray_ = NULL;
866    delete solver;
867    return returnCode;
868}
869// See if diving will give better solution
870// Sets value of solution
871// Returns 1 if solution, 0 if not
872int
873CbcHeuristicDive::solution(double & solutionValue,
874                           double * betterSolution)
875{
876    int nodeCount = model_->getNodeCount();
877    if (feasibilityPumpOptions_>0 && (nodeCount % feasibilityPumpOptions_) != 0)
878        return 0;
879#ifdef DIVE_DEBUG
880    std::cout << "solutionValue = " << solutionValue << std::endl;
881#endif
882    ++numCouldRun_;
883
884    // test if the heuristic can run
885    if (!canHeuristicRun())
886        return 0;
887
888#ifdef JJF_ZERO
889    // See if to do
890    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
891            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
892        return 0; // switched off
893#endif
894    // Get solution array for heuristic solution
895    int numberColumns = model_->solver()->getNumCols();
896    double * newSolution = new double [numberColumns];
897    int numberCuts=0;
898    int numberNodes=-1;
899    CbcSubProblem ** nodes=NULL;
900    int returnCode=solution(solutionValue,numberNodes,numberCuts,
901                            NULL,nodes,
902                            newSolution);
903  if (returnCode==1) 
904    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
905
906    delete [] newSolution;
907    return returnCode;
908}
909/* returns 0 if no solution, 1 if valid solution
910   with better objective value than one passed in
911   also returns list of nodes
912   This does Fractional Diving
913*/
914int 
915CbcHeuristicDive::fathom(CbcModel * model, int & numberNodes, 
916                         CbcSubProblem ** & nodes)
917{
918  double solutionValue = model->getCutoff();
919  numberNodes=0;
920  // Get solution array for heuristic solution
921  int numberColumns = model_->solver()->getNumCols();
922  double * newSolution = new double [4*numberColumns];
923  double * lastDjs = newSolution+numberColumns;
924  double * originalLower = lastDjs+numberColumns;
925  double * originalUpper = originalLower+numberColumns;
926  memcpy(originalLower,model_->solver()->getColLower(),
927         numberColumns*sizeof(double));
928  memcpy(originalUpper,model_->solver()->getColUpper(),
929         numberColumns*sizeof(double));
930  int numberCuts=0;
931  OsiRowCut ** cuts = NULL; //new OsiRowCut * [maxIterations_];
932  nodes=new CbcSubProblem * [maxIterations_+2];
933  int returnCode=solution(solutionValue,numberNodes,numberCuts,
934                          cuts,nodes,
935                          newSolution);
936
937  if (returnCode==1) {
938    // copy to best solution ? or put in solver
939    printf("Solution from heuristic fathom\n");
940  }
941  int numberFeasibleNodes=numberNodes;
942  if (returnCode!=1)
943    numberFeasibleNodes--;
944  if (numberFeasibleNodes>0) {
945    CoinWarmStartBasis * basis = nodes[numberFeasibleNodes-1]->status_;
946    //double * sort = new double [numberFeasibleNodes];
947    //int * whichNode = new int [numberFeasibleNodes];
948    //int numberNodesNew=0;
949    // use djs on previous unless feasible
950    for (int iNode=0;iNode<numberFeasibleNodes;iNode++) {
951      CbcSubProblem * sub = nodes[iNode];
952      double branchValue = sub->branchValue_;
953      int iStatus=sub->problemStatus_;
954      int iColumn = sub->branchVariable_;
955      bool secondBranch = (iStatus&2)!=0;
956      bool branchUp;
957      if (!secondBranch)
958        branchUp = (iStatus&1)!=0;
959      else
960        branchUp = (iStatus&1)==0;
961      double djValue=lastDjs[iColumn];
962      sub->djValue_=fabs(djValue);
963      if (!branchUp&&floor(branchValue)==originalLower[iColumn]
964          &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atLowerBound) {
965        if (djValue>0.0) {
966          // naturally goes to LB
967          printf("ignoring branch down on %d (node %d) from value of %g - branch was %s - dj %g\n",
968                 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
969                 djValue);
970          sub->problemStatus_ |= 4;
971          //} else {
972          // put on list
973          //sort[numberNodesNew]=djValue;
974          //whichNode[numberNodesNew++]=iNode;
975        }
976      } else if (branchUp&&ceil(branchValue)==originalUpper[iColumn]
977          &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atUpperBound) {
978        if (djValue<0.0) {
979          // naturally goes to UB
980          printf("ignoring branch up on %d (node %d) from value of %g - branch was %s - dj %g\n",
981                 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
982                 djValue);
983          sub->problemStatus_ |= 4;
984          //} else {
985          // put on list
986          //sort[numberNodesNew]=-djValue;
987          //whichNode[numberNodesNew++]=iNode;
988        }
989      }
990    }
991    // use conflict to order nodes
992    for (int iCut=0;iCut<numberCuts;iCut++) {
993    }
994    //CoinSort_2(sort,sort+numberNodesNew,whichNode);
995    // create nodes
996    // last node will have one way already done
997  }
998  for (int iCut=0;iCut<numberCuts;iCut++) {
999    delete cuts[iCut];
1000  }
1001  delete [] cuts;
1002  delete [] newSolution;
1003  return returnCode;
1004}
1005
1006// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1007void
1008CbcHeuristicDive::validate()
1009{
1010    if (model_ && (when() % 100) < 10) {
1011        if (model_->numberIntegers() !=
1012                model_->numberObjects() && (model_->numberObjects() ||
1013                                            (model_->specialOptions()&1024) == 0)) {
1014            int numberOdd = 0;
1015            for (int i = 0; i < model_->numberObjects(); i++) {
1016                if (!model_->object(i)->canDoHeuristics())
1017                    numberOdd++;
1018            }
1019            if (numberOdd)
1020                setWhen(0);
1021        }
1022    }
1023
1024    int numberIntegers = model_->numberIntegers();
1025    const int * integerVariable = model_->integerVariable();
1026    delete [] downLocks_;
1027    delete [] upLocks_;
1028    downLocks_ = new unsigned short [numberIntegers];
1029    upLocks_ = new unsigned short [numberIntegers];
1030    // Column copy
1031    const double * element = matrix_.getElements();
1032    const int * row = matrix_.getIndices();
1033    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1034    const int * columnLength = matrix_.getVectorLengths();
1035    const double * rowLower = model_->solver()->getRowLower();
1036    const double * rowUpper = model_->solver()->getRowUpper();
1037    for (int i = 0; i < numberIntegers; i++) {
1038        int iColumn = integerVariable[i];
1039        int down = 0;
1040        int up = 0;
1041        if (columnLength[iColumn] > 65535) {
1042            setWhen(0);
1043            break; // unlikely to work
1044        }
1045        for (CoinBigIndex j = columnStart[iColumn];
1046                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1047            int iRow = row[j];
1048            if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
1049                up++;
1050                down++;
1051            } else if (element[j] > 0.0) {
1052                if (rowUpper[iRow] < 1.0e20)
1053                    up++;
1054                else
1055                    down++;
1056            } else {
1057                if (rowLower[iRow] > -1.0e20)
1058                    up++;
1059                else
1060                    down++;
1061            }
1062        }
1063        downLocks_[i] = static_cast<unsigned short> (down);
1064        upLocks_[i] = static_cast<unsigned short> (up);
1065    }
1066
1067#ifdef DIVE_FIX_BINARY_VARIABLES
1068    selectBinaryVariables();
1069#endif
1070}
1071
1072// Select candidate binary variables for fixing
1073void
1074CbcHeuristicDive::selectBinaryVariables()
1075{
1076    // Row copy
1077    const double * elementByRow = matrixByRow_.getElements();
1078    const int * column = matrixByRow_.getIndices();
1079    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
1080    const int * rowLength = matrixByRow_.getVectorLengths();
1081
1082    const int numberRows = matrixByRow_.getNumRows();
1083    const int numberCols = matrixByRow_.getNumCols();
1084
1085    const double * lower = model_->solver()->getColLower();
1086    const double * upper = model_->solver()->getColUpper();
1087    const double * rowLower = model_->solver()->getRowLower();
1088    const double * rowUpper = model_->solver()->getRowUpper();
1089
1090    //  const char * integerType = model_->integerType();
1091
1092
1093    //  const int numberIntegers = model_->numberIntegers();
1094    //  const int * integerVariable = model_->integerVariable();
1095    const double * objective = model_->solver()->getObjCoefficients();
1096
1097    // vector to store the row number of variable bound rows
1098    int* rowIndexes = new int [numberCols];
1099    memset(rowIndexes, -1, numberCols*sizeof(int));
1100
1101    for (int i = 0; i < numberRows; i++) {
1102        int positiveBinary = -1;
1103        int negativeBinary = -1;
1104        int nPositiveOther = 0;
1105        int nNegativeOther = 0;
1106        for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
1107            int iColumn = column[k];
1108            if (model_->solver()->isInteger(iColumn) &&
1109                    lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
1110                    objective[iColumn] == 0.0 &&
1111                    elementByRow[k] > 0.0 &&
1112                    positiveBinary < 0)
1113                positiveBinary = iColumn;
1114            else if (model_->solver()->isInteger(iColumn) &&
1115                     lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
1116                     objective[iColumn] == 0.0 &&
1117                     elementByRow[k] < 0.0 &&
1118                     negativeBinary < 0)
1119                negativeBinary = iColumn;
1120            else if ((elementByRow[k] > 0.0 &&
1121                      lower[iColumn] >= 0.0) ||
1122                     (elementByRow[k] < 0.0 &&
1123                      upper[iColumn] <= 0.0))
1124                nPositiveOther++;
1125            else if ((elementByRow[k] > 0.0 &&
1126                      lower[iColumn] <= 0.0) ||
1127                     (elementByRow[k] < 0.0 &&
1128                      upper[iColumn] >= 0.0))
1129                nNegativeOther++;
1130            if (nPositiveOther > 0 && nNegativeOther > 0)
1131                break;
1132        }
1133        int binVar = -1;
1134        if (positiveBinary >= 0 &&
1135                (negativeBinary >= 0 || nNegativeOther > 0) &&
1136                nPositiveOther == 0 &&
1137                rowLower[i] == 0.0 &&
1138                rowUpper[i] > 0.0)
1139            binVar = positiveBinary;
1140        else if (negativeBinary >= 0 &&
1141                 (positiveBinary >= 0 || nPositiveOther > 0) &&
1142                 nNegativeOther == 0 &&
1143                 rowLower[i] < 0.0 &&
1144                 rowUpper[i] == 0.0)
1145            binVar = negativeBinary;
1146        if (binVar >= 0) {
1147            if (rowIndexes[binVar] == -1)
1148                rowIndexes[binVar] = i;
1149            else if (rowIndexes[binVar] >= 0)
1150                rowIndexes[binVar] = -2;
1151        }
1152    }
1153
1154    for (int j = 0; j < numberCols; j++) {
1155        if (rowIndexes[j] >= 0) {
1156            binVarIndex_.push_back(j);
1157            vbRowIndex_.push_back(rowIndexes[j]);
1158        }
1159    }
1160
1161#ifdef DIVE_DEBUG
1162    std::cout << "number vub Binary = " << binVarIndex_.size() << std::endl;
1163#endif
1164
1165    delete [] rowIndexes;
1166
1167}
1168
1169/*
1170  Perform reduced cost fixing on integer variables.
1171
1172  The variables in question are already nonbasic at bound. We're just nailing
1173  down the current situation.
1174*/
1175
1176int CbcHeuristicDive::reducedCostFix (OsiSolverInterface* solver)
1177
1178{
1179    //return 0; // temp
1180#ifndef JJF_ONE
1181    if (!model_->solverCharacteristics()->reducedCostsAccurate())
1182        return 0; //NLP
1183#endif
1184    double cutoff = model_->getCutoff() ;
1185    if (cutoff > 1.0e20)
1186        return 0;
1187#ifdef DIVE_DEBUG
1188    std::cout << "cutoff = " << cutoff << std::endl;
1189#endif
1190    double direction = solver->getObjSense() ;
1191    double gap = cutoff - solver->getObjValue() * direction ;
1192    gap *= 0.5; // Fix more
1193    double tolerance;
1194    solver->getDblParam(OsiDualTolerance, tolerance) ;
1195    if (gap <= 0.0)
1196        gap = tolerance; //return 0;
1197    gap += 100.0 * tolerance;
1198    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1199
1200    const double *lower = solver->getColLower() ;
1201    const double *upper = solver->getColUpper() ;
1202    const double *solution = solver->getColSolution() ;
1203    const double *reducedCost = solver->getReducedCost() ;
1204
1205    int numberIntegers = model_->numberIntegers();
1206    const int * integerVariable = model_->integerVariable();
1207
1208    int numberFixed = 0 ;
1209
1210# ifdef COIN_HAS_CLP
1211    OsiClpSolverInterface * clpSolver
1212    = dynamic_cast<OsiClpSolverInterface *> (solver);
1213    ClpSimplex * clpSimplex = NULL;
1214    if (clpSolver)
1215        clpSimplex = clpSolver->getModelPtr();
1216# endif
1217    for (int i = 0 ; i < numberIntegers ; i++) {
1218        int iColumn = integerVariable[i] ;
1219        double djValue = direction * reducedCost[iColumn] ;
1220        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1221            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1222#ifdef COIN_HAS_CLP
1223                // may just have been fixed before
1224                if (clpSimplex) {
1225                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
1226#ifdef COIN_DEVELOP
1227                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
1228                               iColumn, clpSimplex->getColumnStatus(iColumn),
1229                               djValue, gap, lower[iColumn], upper[iColumn]);
1230#endif
1231                    } else {
1232                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atLowerBound ||
1233                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
1234                    }
1235                }
1236#endif
1237                solver->setColUpper(iColumn, lower[iColumn]) ;
1238                numberFixed++ ;
1239            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1240#ifdef COIN_HAS_CLP
1241                // may just have been fixed before
1242                if (clpSimplex) {
1243                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
1244#ifdef COIN_DEVELOP
1245                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
1246                               iColumn, clpSimplex->getColumnStatus(iColumn),
1247                               djValue, gap, lower[iColumn], upper[iColumn]);
1248#endif
1249                    } else {
1250                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atUpperBound ||
1251                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
1252                    }
1253                }
1254#endif
1255                solver->setColLower(iColumn, upper[iColumn]) ;
1256                numberFixed++ ;
1257            }
1258        }
1259    }
1260    return numberFixed;
1261}
1262// Fix other variables at bounds
1263int
1264CbcHeuristicDive::fixOtherVariables(OsiSolverInterface * solver,
1265                                    const double * solution,
1266                                    PseudoReducedCost * candidate,
1267                                    const double * random)
1268{
1269    const double * lower = solver->getColLower();
1270    const double * upper = solver->getColUpper();
1271    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1272    double primalTolerance;
1273    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1274
1275    int numberIntegers = model_->numberIntegers();
1276    const int * integerVariable = model_->integerVariable();
1277    const double* reducedCost = solver->getReducedCost();
1278    // fix other integer variables that are at their bounds
1279    int cnt = 0;
1280#ifdef GAP
1281    double direction = solver->getObjSense(); // 1 for min, -1 for max
1282    double gap = 1.0e30;
1283#endif
1284#ifdef GAP
1285    double cutoff = model_->getCutoff() ;
1286    if (cutoff < 1.0e20 && false) {
1287        double direction = solver->getObjSense() ;
1288        gap = cutoff - solver->getObjValue() * direction ;
1289        gap *= 0.1; // Fix more if plausible
1290        double tolerance;
1291        solver->getDblParam(OsiDualTolerance, tolerance) ;
1292        if (gap <= 0.0)
1293            gap = tolerance;
1294        gap += 100.0 * tolerance;
1295    }
1296    int nOverGap = 0;
1297#endif
1298    int numberFree = 0;
1299    int numberFixedAlready = 0;
1300    for (int i = 0; i < numberIntegers; i++) {
1301        int iColumn = integerVariable[i];
1302        if (upper[iColumn] > lower[iColumn]) {
1303            numberFree++;
1304            double value = solution[iColumn];
1305            if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
1306                candidate[cnt].var = iColumn;
1307                candidate[cnt++].pseudoRedCost =
1308                    fabs(reducedCost[iColumn] * random[i]);
1309#ifdef GAP
1310                if (fabs(reducedCost[iColumn]) > gap)
1311                    nOverGap++;
1312#endif
1313            }
1314        } else {
1315            numberFixedAlready++;
1316        }
1317    }
1318#ifdef GAP
1319    int nLeft = maxNumberToFix - numberFixedAlready;
1320#ifdef CLP_INVESTIGATE4
1321    printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
1322           cutoff, solver->getObjValue(), nOverGap, numberFree,
1323           numberFixedAlready);
1324#endif
1325    if (nOverGap > nLeft && true) {
1326        nOverGap = CoinMin(nOverGap, nLeft + maxNumberToFix / 2);
1327        maxNumberToFix += nOverGap - nLeft;
1328    }
1329#else
1330#ifdef CLP_INVESTIGATE4
1331    printf("cutoff %g obj %g - %d free, %d fixed\n",
1332           model_->getCutoff(), solver->getObjValue(), numberFree,
1333           numberFixedAlready);
1334#endif
1335#endif
1336    return cnt;
1337}
1338
Note: See TracBrowser for help on using the repository browser.