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

Last change on this file since 2344 was 2344, checked in by forrest, 2 years ago

change int to CoinBigIndex?

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