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

Last change on this file since 2092 was 2092, checked in by forrest, 5 years ago

symmetry and diving

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.7 KB
Line 
1/* $Id: CbcHeuristicDive.cpp 2092 2014-10-03 11:58:33Z forrest $ */
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        // clean
294        value = CoinMin(value,upperBefore[iColumn]);
295        value = CoinMax(value,lowerBefore[iColumn]);
296        newSolution[iColumn] = value;
297        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
298            numberFractionalVariables++;
299        }
300    }
301
302    const double* reducedCost = NULL;
303    // See if not NLP
304    if (model_->solverCharacteristics()->reducedCostsAccurate())
305        reducedCost = solver->getReducedCost();
306
307    int iteration = 0;
308    while (numberFractionalVariables) {
309        iteration++;
310
311        // initialize any data
312        initializeData();
313
314        // select a fractional variable to bound
315        int bestColumn = -1;
316        int bestRound; // -1 rounds down, +1 rounds up
317        bool canRound = selectVariableToBranch(solver, newSolution,
318                                               bestColumn, bestRound);
319        // if the solution is not trivially roundable, we don't try to round;
320        // if the solution is trivially roundable, we try to round. However,
321        // if the rounded solution is worse than the current incumbent,
322        // then we don't round and proceed normally. In this case, the
323        // bestColumn will be a trivially roundable variable
324        if (canRound) {
325            // check if by rounding all fractional variables
326            // we get a solution with an objective value
327            // better than the current best integer solution
328            double delta = 0.0;
329            for (int i = 0; i < numberIntegers; i++) {
330                int iColumn = integerVariable[i];
331                double value = newSolution[iColumn];
332                if (fabs(floor(value + 0.5) - value) > integerTolerance) {
333                    assert(downLocks_[i] == 0 || upLocks_[i] == 0);
334                    double obj = objective[iColumn];
335                    if (downLocks_[i] == 0 && upLocks_[i] == 0) {
336                        if (direction * obj >= 0.0)
337                            delta += (floor(value) - value) * obj;
338                        else
339                            delta += (ceil(value) - value) * obj;
340                    } else if (downLocks_[i] == 0)
341                        delta += (floor(value) - value) * obj;
342                    else
343                        delta += (ceil(value) - value) * obj;
344                }
345            }
346            if (direction*(solver->getObjValue() + delta) < solutionValue) {
347#ifdef DIVE_DEBUG
348                nRoundFeasible++;
349#endif
350                if (!nodes||bestColumn<0) {
351                  // Round all the fractional variables
352                  for (int i = 0; i < numberIntegers; i++) {
353                    int iColumn = integerVariable[i];
354                    double value = newSolution[iColumn];
355                    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
356                      assert(downLocks_[i] == 0 || upLocks_[i] == 0);
357                      if (downLocks_[i] == 0 && upLocks_[i] == 0) {
358                        if (direction * objective[iColumn] >= 0.0)
359                          newSolution[iColumn] = floor(value);
360                        else
361                          newSolution[iColumn] = ceil(value);
362                      } else if (downLocks_[i] == 0)
363                        newSolution[iColumn] = floor(value);
364                      else
365                        newSolution[iColumn] = ceil(value);
366                    }
367                  }
368                  break;
369                } else {
370                  // can't round if going to use in branching
371                  int i;
372                  for (i = 0; i < numberIntegers; i++) {
373                    int iColumn = integerVariable[i];
374                    double value = newSolution[bestColumn];
375                    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
376                      if (iColumn==bestColumn) {
377                        assert(downLocks_[i] == 0 || upLocks_[i] == 0);
378                        double obj = objective[bestColumn];
379                        if (downLocks_[i] == 0 && upLocks_[i] == 0) {
380                          if (direction * obj >= 0.0)
381                            bestRound=-1;
382                          else
383                            bestRound=1;
384                        } else if (downLocks_[i] == 0)
385                          bestRound=-1;
386                        else
387                          bestRound=1;
388                        break;
389                      }
390                    }
391                  }
392                }
393            }
394#ifdef DIVE_DEBUG
395            else
396                nRoundInfeasible++;
397#endif
398        }
399
400        // do reduced cost fixing
401#ifdef DIVE_DEBUG
402        int numberFixed = reducedCostFix(solver);
403        std::cout << "numberReducedCostFixed = " << numberFixed << std::endl;
404#else
405        reducedCostFix(solver);
406#endif
407
408        int numberAtBoundFixed = 0;
409#ifdef DIVE_FIX_BINARY_VARIABLES
410        // fix binary variables based on pseudo reduced cost
411        if (binVarIndex_.size()) {
412            int cnt = 0;
413            int n = static_cast<int>(binVarIndex_.size());
414            for (int j = 0; j < n; j++) {
415                int iColumn1 = binVarIndex_[j];
416                double value = newSolution[iColumn1];
417                if (fabs(value) <= integerTolerance &&
418                        lower[iColumn1] != upper[iColumn1]) {
419                    double maxPseudoReducedCost = 0.0;
420#ifdef DIVE_DEBUG
421                    std::cout << "iColumn1 = " << iColumn1 << ", value = " << value << std::endl;
422#endif
423                    int iRow = vbRowIndex_[j];
424                    double chosenValue = 0.0;
425                    for (int k = rowStart[iRow]; k < rowStart[iRow] + rowLength[iRow]; k++) {
426                        int iColumn2 = column[k];
427#ifdef DIVE_DEBUG
428                        std::cout << "iColumn2 = " << iColumn2 << std::endl;
429#endif
430                        if (iColumn1 != iColumn2) {
431                            double pseudoReducedCost = fabs(reducedCost[iColumn2] *
432                                                            elementByRow[k]);
433#ifdef DIVE_DEBUG
434                            int k2;
435                            for (k2 = rowStart[iRow]; k2 < rowStart[iRow] + rowLength[iRow]; k2++) {
436                                if (column[k2] == iColumn1)
437                                    break;
438                            }
439                            std::cout << "reducedCost[" << iColumn2 << "] = "
440                                      << reducedCost[iColumn2]
441                                      << ", elementByRow[" << iColumn2 << "] = " << elementByRow[k]
442                                      << ", elementByRow[" << iColumn1 << "] = " << elementByRow[k2]
443                                      << ", pseudoRedCost = " << pseudoReducedCost
444                                      << std::endl;
445#endif
446                            if (pseudoReducedCost > maxPseudoReducedCost)
447                                maxPseudoReducedCost = pseudoReducedCost;
448                        } else {
449                            // save value
450                            chosenValue = fabs(elementByRow[k]);
451                        }
452                    }
453                    assert (chosenValue);
454                    maxPseudoReducedCost /= chosenValue;
455#ifdef DIVE_DEBUG
456                    std::cout << ", maxPseudoRedCost = " << maxPseudoReducedCost << std::endl;
457#endif
458                    candidate[cnt].var = iColumn1;
459                    candidate[cnt++].pseudoRedCost = maxPseudoReducedCost;
460                }
461            }
462#ifdef DIVE_DEBUG
463            std::cout << "candidates for rounding = " << cnt << std::endl;
464#endif
465            std::sort(candidate, candidate + cnt, compareBinaryVars);
466            for (int i = 0; i < cnt; i++) {
467                int iColumn = candidate[i].var;
468                if (numberAtBoundFixed < maxNumberAtBoundToFix) {
469                    columnFixed[numberAtBoundFixed] = iColumn;
470                    originalBound[numberAtBoundFixed] = upper[iColumn];
471                    fixedAtLowerBound[numberAtBoundFixed] = true;
472                    solver->setColUpper(iColumn, lower[iColumn]);
473                    numberAtBoundFixed++;
474                    if (numberAtBoundFixed == maxNumberAtBoundToFix)
475                        break;
476                }
477            }
478        }
479#endif
480
481        // fix other integer variables that are at their bounds
482        int cnt = 0;
483#ifdef GAP
484        double gap = 1.0e30;
485#endif
486        if (reducedCost && true) {
487#ifndef JJF_ONE
488            cnt = fixOtherVariables(solver, solution, candidate, random);
489#else
490#ifdef GAP
491            double cutoff = model_->getCutoff() ;
492            if (cutoff < 1.0e20 && false) {
493                double direction = solver->getObjSense() ;
494                gap = cutoff - solver->getObjValue() * direction ;
495                gap *= 0.1; // Fix more if plausible
496                double tolerance;
497                solver->getDblParam(OsiDualTolerance, tolerance) ;
498                if (gap <= 0.0)
499                    gap = tolerance;
500                gap += 100.0 * tolerance;
501            }
502            int nOverGap = 0;
503#endif
504            int numberFree = 0;
505            int numberFixed = 0;
506            for (int i = 0; i < numberIntegers; i++) {
507                int iColumn = integerVariable[i];
508                if (upper[iColumn] > lower[iColumn]) {
509                    numberFree++;
510                    double value = newSolution[iColumn];
511                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
512                        candidate[cnt].var = iColumn;
513                        candidate[cnt++].pseudoRedCost =
514                            fabs(reducedCost[iColumn] * random[i]);
515#ifdef GAP
516                        if (fabs(reducedCost[iColumn]) > gap)
517                            nOverGap++;
518#endif
519                    }
520                } else {
521                    numberFixed++;
522                }
523            }
524#ifdef GAP
525            int nLeft = maxNumberAtBoundToFix - numberAtBoundFixed;
526#ifdef CLP_INVESTIGATE4
527            printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
528                   cutoff, solver->getObjValue(), nOverGap, numberFree, numberFixed);
529#endif
530            if (nOverGap > nLeft && true) {
531                nOverGap = CoinMin(nOverGap, nLeft + maxNumberAtBoundToFix / 2);
532                maxNumberAtBoundToFix += nOverGap - nLeft;
533            }
534#else
535#ifdef CLP_INVESTIGATE4
536            printf("cutoff %g obj %g - %d free, %d fixed\n",
537                   model_->getCutoff(), solver->getObjValue(), numberFree, numberFixed);
538#endif
539#endif
540#endif
541        } else {
542            for (int i = 0; i < numberIntegers; i++) {
543                int iColumn = integerVariable[i];
544                if (upper[iColumn] > lower[iColumn]) {
545                    double value = newSolution[iColumn];
546                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
547                        candidate[cnt].var = iColumn;
548                        candidate[cnt++].pseudoRedCost = numberIntegers - i;
549                    }
550                }
551            }
552        }
553        std::sort(candidate, candidate + cnt, compareBinaryVars);
554        for (int i = 0; i < cnt; i++) {
555            int iColumn = candidate[i].var;
556            if (upper[iColumn] > lower[iColumn]) {
557                double value = newSolution[iColumn];
558                if (fabs(floor(value + 0.5) - value) <= integerTolerance &&
559                        numberAtBoundFixed < maxNumberAtBoundToFix) {
560                    // fix the variable at one of its bounds
561                    if (fabs(lower[iColumn] - value) <= integerTolerance) {
562                        columnFixed[numberAtBoundFixed] = iColumn;
563                        originalBound[numberAtBoundFixed] = upper[iColumn];
564                        fixedAtLowerBound[numberAtBoundFixed] = true;
565                        solver->setColUpper(iColumn, lower[iColumn]);
566                        numberAtBoundFixed++;
567                    } else if (fabs(upper[iColumn] - value) <= integerTolerance) {
568                        columnFixed[numberAtBoundFixed] = iColumn;
569                        originalBound[numberAtBoundFixed] = lower[iColumn];
570                        fixedAtLowerBound[numberAtBoundFixed] = false;
571                        solver->setColLower(iColumn, upper[iColumn]);
572                        numberAtBoundFixed++;
573                    }
574                    if (numberAtBoundFixed == maxNumberAtBoundToFix)
575                        break;
576                }
577            }
578        }
579#ifdef DIVE_DEBUG
580        std::cout << "numberAtBoundFixed = " << numberAtBoundFixed << std::endl;
581#endif
582
583        double originalBoundBestColumn;
584        double bestColumnValue;
585        int whichWay;
586        if (bestColumn >= 0) {
587            bestColumnValue = newSolution[bestColumn];
588            if (bestRound < 0) {
589                originalBoundBestColumn = upper[bestColumn];
590                solver->setColUpper(bestColumn, floor(bestColumnValue));
591                whichWay=0;
592            } else {
593                originalBoundBestColumn = lower[bestColumn];
594                solver->setColLower(bestColumn, ceil(bestColumnValue));
595                whichWay=1;
596            }
597        } else {
598            break;
599        }
600        int originalBestRound = bestRound;
601        int saveModelOptions = model_->specialOptions();
602       
603        while (1) {
604
605            model_->setSpecialOptions(saveModelOptions | 2048);
606            solver->resolve();
607            model_->setSpecialOptions(saveModelOptions);
608            if (!solver->isAbandoned()&&!solver->isIterationLimitReached()) {
609                numberSimplexIterations += solver->getIterationCount();
610            } else {
611                numberSimplexIterations = maxSimplexIterations + 1;
612                reasonToStop += 100;
613                break;
614            }
615
616            if (!solver->isProvenOptimal()) {
617                if (nodes) {
618                  if (solver->isProvenPrimalInfeasible()) {
619                    if (maxSimplexIterationsAtRoot_!=COIN_INT_MAX) {
620                      // stop now
621                      printf("stopping on first infeasibility\n");
622                      break;
623                    } else if (cuts) {
624                      // can do conflict cut
625                      printf("could do intermediate conflict cut\n");
626                      bool localCut;
627                      OsiRowCut * cut = model_->conflictCut(solver,localCut);
628                      if (cut) {
629                        if (!localCut) {
630                          model_->makePartialCut(cut,solver);
631                          cuts[numberCuts++]=cut;
632                        } else {
633                          delete cut;
634                        }
635                      }
636                    }
637                  } else {
638                    reasonToStop += 10;
639                    break;
640                  }
641                }
642                if (numberAtBoundFixed > 0) {
643                    // Remove the bound fix for variables that were at bounds
644                    for (int i = 0; i < numberAtBoundFixed; i++) {
645                        int iColFixed = columnFixed[i];
646                        if (fixedAtLowerBound[i])
647                            solver->setColUpper(iColFixed, originalBound[i]);
648                        else
649                            solver->setColLower(iColFixed, originalBound[i]);
650                    }
651                    numberAtBoundFixed = 0;
652                } else if (bestRound == originalBestRound) {
653                    bestRound *= (-1);
654                    whichWay |=2;
655                    if (bestRound < 0) {
656                        solver->setColLower(bestColumn, originalBoundBestColumn);
657                        solver->setColUpper(bestColumn, floor(bestColumnValue));
658                    } else {
659                        solver->setColLower(bestColumn, ceil(bestColumnValue));
660                        solver->setColUpper(bestColumn, originalBoundBestColumn);
661                    }
662                } else
663                    break;
664            } else
665                break;
666        }
667
668        if (!solver->isProvenOptimal() ||
669                direction*solver->getObjValue() >= solutionValue) {
670            reasonToStop += 1;
671        } else if (iteration > maxIterations_) {
672            reasonToStop += 2;
673        } else if (CoinCpuTime() - time1 > maxTime_) {
674            reasonToStop += 3;
675        } else if (numberSimplexIterations > maxSimplexIterations) {
676            reasonToStop += 4;
677            // also switch off
678#ifdef CLP_INVESTIGATE
679            printf("switching off diving as too many iterations %d, %d allowed\n",
680                   numberSimplexIterations, maxSimplexIterations);
681#endif
682            when_ = 0;
683        } else if (solver->getIterationCount() > 1000 && iteration > 3 && !nodes) {
684            reasonToStop += 5;
685            // also switch off
686#ifdef CLP_INVESTIGATE
687            printf("switching off diving one iteration took %d iterations (total %d)\n",
688                   solver->getIterationCount(), numberSimplexIterations);
689#endif
690            when_ = 0;
691        }
692
693        memcpy(newSolution, solution, numberColumns*sizeof(double));
694        numberFractionalVariables = 0;
695        double sumFractionalVariables=0.0;
696        for (int i = 0; i < numberIntegers; i++) {
697            int iColumn = integerVariable[i];
698            double value = newSolution[iColumn];
699            double away = fabs(floor(value + 0.5) - value);
700            if (away > integerTolerance) {
701                numberFractionalVariables++;
702                sumFractionalVariables += away;
703            }
704        }
705        if (nodes) {
706          // save information
707          //branchValues[numberNodes]=bestColumnValue;
708          //statuses[numberNodes]=whichWay+(bestColumn<<2);
709          //bases[numberNodes]=solver->getWarmStart();
710          ClpSimplex * simplex = clpSolver->getModelPtr();
711          CbcSubProblem * sub =
712            new CbcSubProblem(clpSolver,lowerBefore,upperBefore,
713                          simplex->statusArray(),numberNodes);
714          nodes[numberNodes]=sub;
715          // other stuff
716          sub->branchValue_=bestColumnValue;
717          sub->problemStatus_=whichWay;
718          sub->branchVariable_=bestColumn;
719          sub->objectiveValue_ = simplex->objectiveValue();
720          sub->sumInfeasibilities_ = sumFractionalVariables;
721          sub->numberInfeasibilities_ = numberFractionalVariables;
722          printf("DiveNode %d column %d way %d bvalue %g obj %g\n",
723                 numberNodes,sub->branchVariable_,sub->problemStatus_,
724                 sub->branchValue_,sub->objectiveValue_);
725          numberNodes++;
726          if (solver->isProvenOptimal()) {
727            memcpy(lastDjs,solver->getReducedCost(),numberColumns*sizeof(double));
728            memcpy(lowerBefore,lower,numberColumns*sizeof(double));
729            memcpy(upperBefore,upper,numberColumns*sizeof(double));
730          }
731        }
732        if (!numberFractionalVariables||reasonToStop)
733          break;
734    }
735    if (nodes) {
736      printf("Exiting dive for reason %d\n",reasonToStop);
737      if (reasonToStop>1) {
738        printf("problems in diving\n");
739        int whichWay=nodes[numberNodes-1]->problemStatus_;
740        CbcSubProblem * sub;
741        if ((whichWay&2)==0) {
742          // leave both ways
743          sub = new CbcSubProblem(*nodes[numberNodes-1]);
744          nodes[numberNodes++]=sub;
745        } else {
746          sub = nodes[numberNodes-1];
747        }
748        if ((whichWay&1)==0)
749          sub->problemStatus_=whichWay|1;
750        else
751          sub->problemStatus_=whichWay&~1;
752      }
753      if (!numberNodes) {
754        // was good at start! - create fake
755        clpSolver->resolve();
756        ClpSimplex * simplex = clpSolver->getModelPtr();
757        CbcSubProblem * sub =
758          new CbcSubProblem(clpSolver,lowerBefore,upperBefore,
759                            simplex->statusArray(),numberNodes);
760        nodes[numberNodes]=sub;
761        // other stuff
762        sub->branchValue_=0.0;
763        sub->problemStatus_=0;
764        sub->branchVariable_=-1;
765        sub->objectiveValue_ = simplex->objectiveValue();
766        sub->sumInfeasibilities_ = 0.0;
767        sub->numberInfeasibilities_ = 0;
768        printf("DiveNode %d column %d way %d bvalue %g obj %g\n",
769               numberNodes,sub->branchVariable_,sub->problemStatus_,
770               sub->branchValue_,sub->objectiveValue_);
771        numberNodes++;
772        assert (solver->isProvenOptimal());
773      }
774      nodes[numberNodes-1]->problemStatus_ |= 256*reasonToStop;
775      // use djs as well
776      if (solver->isProvenPrimalInfeasible()&&cuts) {
777        // can do conflict cut and re-order
778        printf("could do final conflict cut\n");
779        bool localCut;
780        OsiRowCut * cut = model_->conflictCut(solver,localCut);
781        if (cut) {
782          printf("cut - need to use conflict and previous djs\n");
783          if (!localCut) {
784            model_->makePartialCut(cut,solver);
785            cuts[numberCuts++]=cut;
786          } else {
787            delete cut;
788          }
789        } else {
790          printf("bad conflict - just use previous djs\n");
791        }
792      }
793    }
794   
795    // re-compute new solution value
796    double objOffset = 0.0;
797    solver->getDblParam(OsiObjOffset, objOffset);
798    newSolutionValue = -objOffset;
799    for (int i = 0 ; i < numberColumns ; i++ )
800      newSolutionValue += objective[i] * newSolution[i];
801    newSolutionValue *= direction;
802    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
803    if (newSolutionValue < solutionValue && !reasonToStop) {
804      double * rowActivity = new double[numberRows];
805      memset(rowActivity, 0, numberRows*sizeof(double));
806      // paranoid check
807      memset(rowActivity, 0, numberRows*sizeof(double));
808      for (int i = 0; i < numberColumns; i++) {
809        int j;
810        double value = newSolution[i];
811        if (value) {
812          for (j = columnStart[i];
813               j < columnStart[i] + columnLength[i]; j++) {
814            int iRow = row[j];
815            rowActivity[iRow] += value * element[j];
816          }
817        }
818      }
819      // check was approximately feasible
820      bool feasible = true;
821      for (int i = 0; i < numberRows; i++) {
822        if (rowActivity[i] < rowLower[i]) {
823          if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
824            feasible = false;
825        } else if (rowActivity[i] > rowUpper[i]) {
826          if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
827            feasible = false;
828        }
829      }
830      for (int i = 0; i < numberIntegers; i++) {
831        int iColumn = integerVariable[i];
832        double value = newSolution[iColumn];
833        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
834          feasible = false;
835          break;
836        }
837      }
838      if (feasible) {
839        // new solution
840        solutionValue = newSolutionValue;
841        //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
842        //if (cuts)
843        //clpSolver->getModelPtr()->writeMps("good8.mps", 2);
844        returnCode = 1;
845      } else {
846        // Can easily happen
847        //printf("Debug CbcHeuristicDive giving bad solution\n");
848      }
849      delete [] rowActivity;
850    }
851
852#ifdef DIVE_DEBUG
853    std::cout << "nRoundInfeasible = " << nRoundInfeasible
854              << ", nRoundFeasible = " << nRoundFeasible
855              << ", returnCode = " << returnCode
856              << ", reasonToStop = " << reasonToStop
857              << ", simplexIts = " << numberSimplexIterations
858              << ", iterations = " << iteration << std::endl;
859#endif
860
861    delete [] columnFixed;
862    delete [] originalBound;
863    delete [] fixedAtLowerBound;
864    delete [] candidate;
865    delete [] random;
866    delete [] downArray_;
867    downArray_ = NULL;
868    delete [] upArray_;
869    upArray_ = NULL;
870    delete solver;
871    return returnCode;
872}
873// See if diving will give better solution
874// Sets value of solution
875// Returns 1 if solution, 0 if not
876int
877CbcHeuristicDive::solution(double & solutionValue,
878                           double * betterSolution)
879{
880    int nodeCount = model_->getNodeCount();
881    if (feasibilityPumpOptions_>0 && (nodeCount % feasibilityPumpOptions_) != 0)
882        return 0;
883#ifdef DIVE_DEBUG
884    std::cout << "solutionValue = " << solutionValue << std::endl;
885#endif
886    ++numCouldRun_;
887
888    // test if the heuristic can run
889    if (!canHeuristicRun())
890        return 0;
891
892#ifdef JJF_ZERO
893    // See if to do
894    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
895            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
896        return 0; // switched off
897#endif
898    // Get solution array for heuristic solution
899    int numberColumns = model_->solver()->getNumCols();
900    double * newSolution = new double [numberColumns];
901    int numberCuts=0;
902    int numberNodes=-1;
903    CbcSubProblem ** nodes=NULL;
904    int returnCode=solution(solutionValue,numberNodes,numberCuts,
905                            NULL,nodes,
906                            newSolution);
907  if (returnCode==1) 
908    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
909
910    delete [] newSolution;
911    return returnCode;
912}
913/* returns 0 if no solution, 1 if valid solution
914   with better objective value than one passed in
915   also returns list of nodes
916   This does Fractional Diving
917*/
918int 
919CbcHeuristicDive::fathom(CbcModel * model, int & numberNodes, 
920                         CbcSubProblem ** & nodes)
921{
922  double solutionValue = model->getCutoff();
923  numberNodes=0;
924  // Get solution array for heuristic solution
925  int numberColumns = model_->solver()->getNumCols();
926  double * newSolution = new double [4*numberColumns];
927  double * lastDjs = newSolution+numberColumns;
928  double * originalLower = lastDjs+numberColumns;
929  double * originalUpper = originalLower+numberColumns;
930  memcpy(originalLower,model_->solver()->getColLower(),
931         numberColumns*sizeof(double));
932  memcpy(originalUpper,model_->solver()->getColUpper(),
933         numberColumns*sizeof(double));
934  int numberCuts=0;
935  OsiRowCut ** cuts = NULL; //new OsiRowCut * [maxIterations_];
936  nodes=new CbcSubProblem * [maxIterations_+2];
937  int returnCode=solution(solutionValue,numberNodes,numberCuts,
938                          cuts,nodes,
939                          newSolution);
940
941  if (returnCode==1) {
942    // copy to best solution ? or put in solver
943    printf("Solution from heuristic fathom\n");
944  }
945  int numberFeasibleNodes=numberNodes;
946  if (returnCode!=1)
947    numberFeasibleNodes--;
948  if (numberFeasibleNodes>0) {
949    CoinWarmStartBasis * basis = nodes[numberFeasibleNodes-1]->status_;
950    //double * sort = new double [numberFeasibleNodes];
951    //int * whichNode = new int [numberFeasibleNodes];
952    //int numberNodesNew=0;
953    // use djs on previous unless feasible
954    for (int iNode=0;iNode<numberFeasibleNodes;iNode++) {
955      CbcSubProblem * sub = nodes[iNode];
956      double branchValue = sub->branchValue_;
957      int iStatus=sub->problemStatus_;
958      int iColumn = sub->branchVariable_;
959      bool secondBranch = (iStatus&2)!=0;
960      bool branchUp;
961      if (!secondBranch)
962        branchUp = (iStatus&1)!=0;
963      else
964        branchUp = (iStatus&1)==0;
965      double djValue=lastDjs[iColumn];
966      sub->djValue_=fabs(djValue);
967      if (!branchUp&&floor(branchValue)==originalLower[iColumn]
968          &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atLowerBound) {
969        if (djValue>0.0) {
970          // naturally goes to LB
971          printf("ignoring branch down on %d (node %d) from value of %g - branch was %s - dj %g\n",
972                 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
973                 djValue);
974          sub->problemStatus_ |= 4;
975          //} else {
976          // put on list
977          //sort[numberNodesNew]=djValue;
978          //whichNode[numberNodesNew++]=iNode;
979        }
980      } else if (branchUp&&ceil(branchValue)==originalUpper[iColumn]
981          &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atUpperBound) {
982        if (djValue<0.0) {
983          // naturally goes to UB
984          printf("ignoring branch up on %d (node %d) from value of %g - branch was %s - dj %g\n",
985                 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
986                 djValue);
987          sub->problemStatus_ |= 4;
988          //} else {
989          // put on list
990          //sort[numberNodesNew]=-djValue;
991          //whichNode[numberNodesNew++]=iNode;
992        }
993      }
994    }
995    // use conflict to order nodes
996    for (int iCut=0;iCut<numberCuts;iCut++) {
997    }
998    //CoinSort_2(sort,sort+numberNodesNew,whichNode);
999    // create nodes
1000    // last node will have one way already done
1001  }
1002  for (int iCut=0;iCut<numberCuts;iCut++) {
1003    delete cuts[iCut];
1004  }
1005  delete [] cuts;
1006  delete [] newSolution;
1007  return returnCode;
1008}
1009
1010// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1011void
1012CbcHeuristicDive::validate()
1013{
1014    if (model_ && (when() % 100) < 10) {
1015        if (model_->numberIntegers() !=
1016                model_->numberObjects() && (model_->numberObjects() ||
1017                                            (model_->specialOptions()&1024) == 0)) {
1018            int numberOdd = 0;
1019            for (int i = 0; i < model_->numberObjects(); i++) {
1020                if (!model_->object(i)->canDoHeuristics())
1021                    numberOdd++;
1022            }
1023            if (numberOdd)
1024                setWhen(0);
1025        }
1026    }
1027
1028    int numberIntegers = model_->numberIntegers();
1029    const int * integerVariable = model_->integerVariable();
1030    delete [] downLocks_;
1031    delete [] upLocks_;
1032    downLocks_ = new unsigned short [numberIntegers];
1033    upLocks_ = new unsigned short [numberIntegers];
1034    // Column copy
1035    const double * element = matrix_.getElements();
1036    const int * row = matrix_.getIndices();
1037    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1038    const int * columnLength = matrix_.getVectorLengths();
1039    const double * rowLower = model_->solver()->getRowLower();
1040    const double * rowUpper = model_->solver()->getRowUpper();
1041    for (int i = 0; i < numberIntegers; i++) {
1042        int iColumn = integerVariable[i];
1043        int down = 0;
1044        int up = 0;
1045        if (columnLength[iColumn] > 65535) {
1046            setWhen(0);
1047            break; // unlikely to work
1048        }
1049        for (CoinBigIndex j = columnStart[iColumn];
1050                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1051            int iRow = row[j];
1052            if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
1053                up++;
1054                down++;
1055            } else if (element[j] > 0.0) {
1056                if (rowUpper[iRow] < 1.0e20)
1057                    up++;
1058                else
1059                    down++;
1060            } else {
1061                if (rowLower[iRow] > -1.0e20)
1062                    up++;
1063                else
1064                    down++;
1065            }
1066        }
1067        downLocks_[i] = static_cast<unsigned short> (down);
1068        upLocks_[i] = static_cast<unsigned short> (up);
1069    }
1070
1071#ifdef DIVE_FIX_BINARY_VARIABLES
1072    selectBinaryVariables();
1073#endif
1074}
1075
1076// Select candidate binary variables for fixing
1077void
1078CbcHeuristicDive::selectBinaryVariables()
1079{
1080    // Row copy
1081    const double * elementByRow = matrixByRow_.getElements();
1082    const int * column = matrixByRow_.getIndices();
1083    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
1084    const int * rowLength = matrixByRow_.getVectorLengths();
1085
1086    const int numberRows = matrixByRow_.getNumRows();
1087    const int numberCols = matrixByRow_.getNumCols();
1088
1089    const double * lower = model_->solver()->getColLower();
1090    const double * upper = model_->solver()->getColUpper();
1091    const double * rowLower = model_->solver()->getRowLower();
1092    const double * rowUpper = model_->solver()->getRowUpper();
1093
1094    //  const char * integerType = model_->integerType();
1095
1096
1097    //  const int numberIntegers = model_->numberIntegers();
1098    //  const int * integerVariable = model_->integerVariable();
1099    const double * objective = model_->solver()->getObjCoefficients();
1100
1101    // vector to store the row number of variable bound rows
1102    int* rowIndexes = new int [numberCols];
1103    memset(rowIndexes, -1, numberCols*sizeof(int));
1104
1105    for (int i = 0; i < numberRows; i++) {
1106        int positiveBinary = -1;
1107        int negativeBinary = -1;
1108        int nPositiveOther = 0;
1109        int nNegativeOther = 0;
1110        for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
1111            int iColumn = column[k];
1112            if (model_->solver()->isInteger(iColumn) &&
1113                    lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
1114                    objective[iColumn] == 0.0 &&
1115                    elementByRow[k] > 0.0 &&
1116                    positiveBinary < 0)
1117                positiveBinary = iColumn;
1118            else if (model_->solver()->isInteger(iColumn) &&
1119                     lower[iColumn] == 0.0 && upper[iColumn] == 1.0 &&
1120                     objective[iColumn] == 0.0 &&
1121                     elementByRow[k] < 0.0 &&
1122                     negativeBinary < 0)
1123                negativeBinary = iColumn;
1124            else if ((elementByRow[k] > 0.0 &&
1125                      lower[iColumn] >= 0.0) ||
1126                     (elementByRow[k] < 0.0 &&
1127                      upper[iColumn] <= 0.0))
1128                nPositiveOther++;
1129            else if ((elementByRow[k] > 0.0 &&
1130                      lower[iColumn] <= 0.0) ||
1131                     (elementByRow[k] < 0.0 &&
1132                      upper[iColumn] >= 0.0))
1133                nNegativeOther++;
1134            if (nPositiveOther > 0 && nNegativeOther > 0)
1135                break;
1136        }
1137        int binVar = -1;
1138        if (positiveBinary >= 0 &&
1139                (negativeBinary >= 0 || nNegativeOther > 0) &&
1140                nPositiveOther == 0 &&
1141                rowLower[i] == 0.0 &&
1142                rowUpper[i] > 0.0)
1143            binVar = positiveBinary;
1144        else if (negativeBinary >= 0 &&
1145                 (positiveBinary >= 0 || nPositiveOther > 0) &&
1146                 nNegativeOther == 0 &&
1147                 rowLower[i] < 0.0 &&
1148                 rowUpper[i] == 0.0)
1149            binVar = negativeBinary;
1150        if (binVar >= 0) {
1151            if (rowIndexes[binVar] == -1)
1152                rowIndexes[binVar] = i;
1153            else if (rowIndexes[binVar] >= 0)
1154                rowIndexes[binVar] = -2;
1155        }
1156    }
1157
1158    for (int j = 0; j < numberCols; j++) {
1159        if (rowIndexes[j] >= 0) {
1160            binVarIndex_.push_back(j);
1161            vbRowIndex_.push_back(rowIndexes[j]);
1162        }
1163    }
1164
1165#ifdef DIVE_DEBUG
1166    std::cout << "number vub Binary = " << binVarIndex_.size() << std::endl;
1167#endif
1168
1169    delete [] rowIndexes;
1170
1171}
1172
1173/*
1174  Perform reduced cost fixing on integer variables.
1175
1176  The variables in question are already nonbasic at bound. We're just nailing
1177  down the current situation.
1178*/
1179
1180int CbcHeuristicDive::reducedCostFix (OsiSolverInterface* solver)
1181
1182{
1183    //return 0; // temp
1184#ifndef JJF_ONE
1185    if (!model_->solverCharacteristics()->reducedCostsAccurate())
1186        return 0; //NLP
1187#endif
1188    double cutoff = model_->getCutoff() ;
1189    if (cutoff > 1.0e20)
1190        return 0;
1191#ifdef DIVE_DEBUG
1192    std::cout << "cutoff = " << cutoff << std::endl;
1193#endif
1194    double direction = solver->getObjSense() ;
1195    double gap = cutoff - solver->getObjValue() * direction ;
1196    gap *= 0.5; // Fix more
1197    double tolerance;
1198    solver->getDblParam(OsiDualTolerance, tolerance) ;
1199    if (gap <= 0.0)
1200        gap = tolerance; //return 0;
1201    gap += 100.0 * tolerance;
1202    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1203
1204    const double *lower = solver->getColLower() ;
1205    const double *upper = solver->getColUpper() ;
1206    const double *solution = solver->getColSolution() ;
1207    const double *reducedCost = solver->getReducedCost() ;
1208
1209    int numberIntegers = model_->numberIntegers();
1210    const int * integerVariable = model_->integerVariable();
1211
1212    int numberFixed = 0 ;
1213
1214# ifdef COIN_HAS_CLP
1215    OsiClpSolverInterface * clpSolver
1216    = dynamic_cast<OsiClpSolverInterface *> (solver);
1217    ClpSimplex * clpSimplex = NULL;
1218    if (clpSolver)
1219        clpSimplex = clpSolver->getModelPtr();
1220# endif
1221    for (int i = 0 ; i < numberIntegers ; i++) {
1222        int iColumn = integerVariable[i] ;
1223        double djValue = direction * reducedCost[iColumn] ;
1224        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1225            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1226#ifdef COIN_HAS_CLP
1227                // may just have been fixed before
1228                if (clpSimplex) {
1229                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
1230#ifdef COIN_DEVELOP
1231                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
1232                               iColumn, clpSimplex->getColumnStatus(iColumn),
1233                               djValue, gap, lower[iColumn], upper[iColumn]);
1234#endif
1235                    } else {
1236                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atLowerBound ||
1237                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
1238                    }
1239                }
1240#endif
1241                solver->setColUpper(iColumn, lower[iColumn]) ;
1242                numberFixed++ ;
1243            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1244#ifdef COIN_HAS_CLP
1245                // may just have been fixed before
1246                if (clpSimplex) {
1247                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
1248#ifdef COIN_DEVELOP
1249                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
1250                               iColumn, clpSimplex->getColumnStatus(iColumn),
1251                               djValue, gap, lower[iColumn], upper[iColumn]);
1252#endif
1253                    } else {
1254                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atUpperBound ||
1255                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
1256                    }
1257                }
1258#endif
1259                solver->setColLower(iColumn, upper[iColumn]) ;
1260                numberFixed++ ;
1261            }
1262        }
1263    }
1264    return numberFixed;
1265}
1266// Fix other variables at bounds
1267int
1268CbcHeuristicDive::fixOtherVariables(OsiSolverInterface * solver,
1269                                    const double * solution,
1270                                    PseudoReducedCost * candidate,
1271                                    const double * random)
1272{
1273    const double * lower = solver->getColLower();
1274    const double * upper = solver->getColUpper();
1275    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1276    double primalTolerance;
1277    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1278
1279    int numberIntegers = model_->numberIntegers();
1280    const int * integerVariable = model_->integerVariable();
1281    const double* reducedCost = solver->getReducedCost();
1282    // fix other integer variables that are at their bounds
1283    int cnt = 0;
1284#ifdef GAP
1285    double direction = solver->getObjSense(); // 1 for min, -1 for max
1286    double gap = 1.0e30;
1287#endif
1288#ifdef GAP
1289    double cutoff = model_->getCutoff() ;
1290    if (cutoff < 1.0e20 && false) {
1291        double direction = solver->getObjSense() ;
1292        gap = cutoff - solver->getObjValue() * direction ;
1293        gap *= 0.1; // Fix more if plausible
1294        double tolerance;
1295        solver->getDblParam(OsiDualTolerance, tolerance) ;
1296        if (gap <= 0.0)
1297            gap = tolerance;
1298        gap += 100.0 * tolerance;
1299    }
1300    int nOverGap = 0;
1301#endif
1302    int numberFree = 0;
1303    int numberFixedAlready = 0;
1304    for (int i = 0; i < numberIntegers; i++) {
1305        int iColumn = integerVariable[i];
1306        if (upper[iColumn] > lower[iColumn]) {
1307            numberFree++;
1308            double value = solution[iColumn];
1309            if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
1310                candidate[cnt].var = iColumn;
1311                candidate[cnt++].pseudoRedCost =
1312                    fabs(reducedCost[iColumn] * random[i]);
1313#ifdef GAP
1314                if (fabs(reducedCost[iColumn]) > gap)
1315                    nOverGap++;
1316#endif
1317            }
1318        } else {
1319            numberFixedAlready++;
1320        }
1321    }
1322#ifdef GAP
1323    int nLeft = maxNumberToFix - numberFixedAlready;
1324#ifdef CLP_INVESTIGATE4
1325    printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
1326           cutoff, solver->getObjValue(), nOverGap, numberFree,
1327           numberFixedAlready);
1328#endif
1329    if (nOverGap > nLeft && true) {
1330        nOverGap = CoinMin(nOverGap, nLeft + maxNumberToFix / 2);
1331        maxNumberToFix += nOverGap - nLeft;
1332    }
1333#else
1334#ifdef CLP_INVESTIGATE4
1335    printf("cutoff %g obj %g - %d free, %d fixed\n",
1336           model_->getCutoff(), solver->getObjValue(), numberFree,
1337           numberFixedAlready);
1338#endif
1339#endif
1340    return cnt;
1341}
1342
Note: See TracBrowser for help on using the repository browser.