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

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

changes for diving heuristic

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