source: trunk/Cbc/src/CbcHeuristicDive.cpp

Last change on this file was 2467, checked in by unxusr, 11 months ago

spaces after angles

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