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

Last change on this file since 2187 was 2187, checked in by stefan, 3 years ago

more sync with trunk

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