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

Last change on this file since 2101 was 2101, checked in by forrest, 3 years ago

minor print change to diving

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