source: stable/2.9/Cbc/src/CbcHeuristicDive.cpp @ 2143

Last change on this file since 2143 was 2127, checked in by forrest, 4 years ago

fix for sos before integers in dive

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