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

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

for memory leaks and heuristics and some experimental stuff

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