source: trunk/Cbc/src/CbcHeuristicRENS.cpp @ 1573

Last change on this file since 1573 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

File size: 16.6 KB
RevLine 
[1573]1// $Id$
2// Copyright (C) 2006, 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// edwin 12/5/09 carved out of CbcHeuristicRINS
7
[1368]8#if defined(_MSC_VER)
9// Turn off compiler warning about long names
10#  pragma warning(disable:4786)
11#endif
12#include <cassert>
13#include <cstdlib>
14#include <cmath>
15#include <cfloat>
16
17#include "OsiSolverInterface.hpp"
18#include "CbcModel.hpp"
19#include "CbcMessage.hpp"
20#include "CbcHeuristicRENS.hpp"
[1499]21#include "CoinWarmStartBasis.hpp"
[1564]22#include "CoinSort.hpp"
[1368]23#include "CbcBranchActual.hpp"
24#include "CbcStrategy.hpp"
25#include "CglPreProcess.hpp"
26
27// Default Constructor
28CbcHeuristicRENS::CbcHeuristicRENS()
29        : CbcHeuristic()
30{
31    numberTries_ = 0;
[1499]32    rensType_ = 0;
[1368]33    whereFrom_ = 256 + 1;
34}
35
36// Constructor with model - assumed before cuts
37
38CbcHeuristicRENS::CbcHeuristicRENS(CbcModel & model)
39        : CbcHeuristic(model)
40{
41    numberTries_ = 0;
[1499]42    rensType_ = 0;
[1368]43    whereFrom_ = 256 + 1;
44}
45
46// Destructor
47CbcHeuristicRENS::~CbcHeuristicRENS ()
48{
49}
50
51// Clone
52CbcHeuristic *
53CbcHeuristicRENS::clone() const
54{
55    return new CbcHeuristicRENS(*this);
56}
57
58// Assignment operator
59CbcHeuristicRENS &
60CbcHeuristicRENS::operator=( const CbcHeuristicRENS & rhs)
61{
62    if (this != &rhs) {
63        CbcHeuristic::operator=(rhs);
64        numberTries_ = rhs.numberTries_;
[1499]65        rensType_ = rhs.rensType_;
[1368]66    }
67    return *this;
68}
69
70// Copy constructor
71CbcHeuristicRENS::CbcHeuristicRENS(const CbcHeuristicRENS & rhs)
72        :
73        CbcHeuristic(rhs),
[1499]74        numberTries_(rhs.numberTries_),
75        rensType_(rhs.rensType_)
[1368]76{
77}
78// Resets stuff if model changes
79void
80CbcHeuristicRENS::resetModel(CbcModel * )
81{
82}
83int
84CbcHeuristicRENS::solution(double & solutionValue,
85                           double * betterSolution)
86{
87    int returnCode = 0;
88    const double * bestSolution = model_->bestSolution();
[1499]89    if ((numberTries_&&(rensType_&16)==0) || numberTries_>1 || (when() < 2 && bestSolution))
[1368]90        return 0;
91    numberTries_++;
92    OsiSolverInterface * solver = model_->solver();
93
94    int numberIntegers = model_->numberIntegers();
95    const int * integerVariable = model_->integerVariable();
96
97    const double * currentSolution = solver->getColSolution();
[1569]98    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
99    newSolver->resolve();
100    double direction = newSolver->getObjSense();
101    double cutoff=model_->getCutoff();
102    //newSolver->getDblParam(OsiDualObjectiveLimit, cutoff);
103    //cutoff *= direction;
104    double gap = cutoff - newSolver->getObjValue() * direction ;
105    double tolerance;
106    newSolver->getDblParam(OsiDualTolerance, tolerance) ;
107    if (gap > 0.0 || !newSolver->isProvenOptimal()) {
108      gap += 100.0 * tolerance;
109      int nFix = newSolver->reducedCostFix(gap);
110      if (nFix) {
111        char line [200];
112        sprintf(line, "Reduced cost fixing fixed %d variables", nFix);
113        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
114          << line
115          << CoinMessageEol;
116      }
117    } else {
118      return 0; // finished?
119    }
[1564]120    int numberColumns = solver->getNumCols();
121    double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
[1499]122    int type = rensType_&15;
123    double djTolerance = (type!=1) ? -1.0e30 : 1.0e-4;
[1368]124    const double * colLower = newSolver->getColLower();
125    const double * colUpper = newSolver->getColUpper();
[1564]126    int numberFixed = 0;
[1569]127    if (type==3) {
[1564]128      double total=0.0;
129      int n=0;
130      CoinWarmStartBasis * basis =
131        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
132      if (basis&&basis->getNumArtificial()) {
133        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
134          if (colUpper[iColumn]>colLower[iColumn]&&
135              basis->getStructStatus(iColumn) !=
136              CoinWarmStartBasis::basic) {
137            n++;
138            total += fabs(dj[iColumn]);
139          }
140        }
141        if (n)
142          djTolerance = (0.01*total)/static_cast<double>(n);
143        delete basis;
144      }
[1569]145    } else if (type>=5&&type<=10) {
146      /* 5 fix sets at one
147         6 fix on dj but leave unfixed SOS slacks
148         7 fix sets at one but use pi
149         8 fix all at zero but leave unfixed SOS slacks
150         9 as 8 but only fix all at zero if just one in set nonzero
151         10 as 7 but pi other way
152      */
[1564]153      // SOS type fixing
[1569]154      bool fixSets = (type==5)||(type==7)||(type==10);
[1564]155      CoinWarmStartBasis * basis =
156        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
157      if (basis&&basis->getNumArtificial()) {
158        //const double * rowLower = solver->getRowLower();
159        const double * rowUpper = solver->getRowUpper();
160       
161        int numberRows = solver->getNumRows();
162        // Column copy
163        const CoinPackedMatrix * matrix = solver->getMatrixByCol();
164        const double * element = matrix->getElements();
165        const int * row = matrix->getIndices();
166        const CoinBigIndex * columnStart = matrix->getVectorStarts();
167        const int * columnLength = matrix->getVectorLengths();
168        double * bestDj = new double [numberRows];
169        for (int i=0;i<numberRows;i++) {
170          if (rowUpper[i]==1.0)
171            bestDj[i]=1.0e20;
172          else
173            bestDj[i]=1.0e30;
174        }
175        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
176          if (colUpper[iColumn]>colLower[iColumn]) {
177            CoinBigIndex j;
178            if (currentSolution[iColumn]>1.0e-6&&
179                currentSolution[iColumn]<0.999999) {
180              for (j = columnStart[iColumn];
181                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
182                int iRow = row[j];
[1569]183                if (bestDj[iRow]<1.0e30) {
184                  if (element[j] != 1.0)
185                    bestDj[iRow]=1.0e30;
186                  else
187                    bestDj[iRow]=1.0e25;
188                }
[1499]189              }
[1564]190            } else if ( basis->getStructStatus(iColumn) !=
191              CoinWarmStartBasis::basic) {
192              for (j = columnStart[iColumn];
193                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
194                int iRow = row[j];
[1569]195                if (bestDj[iRow]<1.0e25) {
[1564]196                  if (element[j] != 1.0)
197                    bestDj[iRow]=1.0e30;
198                  else
199                    bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]);
200                }
201              }
[1499]202            }
[1564]203          }
[1499]204        }
[1564]205        int nSOS=0;
206        double * sort = new double [numberRows];
[1569]207        const double * pi = newSolver->getRowPrice();
[1564]208        for (int i=0;i<numberRows;i++) {
209          if (bestDj[i]<1.0e30) {
[1569]210            if (type==5)
211              sort[nSOS++]=bestDj[i];
212            else if (type==7)
213              sort[nSOS++]=-fabs(pi[i]);
214            else
215              sort[nSOS++]=fabs(pi[i]);
[1564]216          }
217        }
218        if (10*nSOS>8*numberRows) {
219          std::sort(sort,sort+nSOS);
220          int last = static_cast<int>(nSOS*0.9*fractionSmall_);
221          double tolerance = sort[last];
222          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
223            if (colUpper[iColumn]>colLower[iColumn]) {
224              CoinBigIndex j;
225              if (currentSolution[iColumn]<=1.0e-6||
226                  currentSolution[iColumn]>=0.999999) {
227                if (fixSets) {
228                  for (j = columnStart[iColumn];
229                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
230                    int iRow = row[j];
[1569]231                    double useDj;
232                    if (type==5) 
233                      useDj = bestDj[iRow];
234                    else if (type==7)
235                      useDj= -fabs(pi[iRow]);
236                    else
237                      useDj= fabs(pi[iRow]);
238                    if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
[1564]239                      numberFixed++;
240                      if (currentSolution[iColumn]<=1.0e-6)
241                        newSolver->setColUpper(iColumn,0.0);
242                      else if (currentSolution[iColumn]>=0.999999) 
243                        newSolver->setColLower(iColumn,1.0);
244                    }
245                  }
246                } else if (columnLength[iColumn]==1) {
247                  // leave more slacks
248                  int iRow = row[columnStart[iColumn]];
249                  if (bestDj[iRow]<1.0e30) {
250                    // fake dj
251                    dj[iColumn] *= 0.000001;
252                  }
[1569]253                } else if (type==8||type==9) {
254                  if (currentSolution[iColumn]<=1.0e-6) {
255                    if (type==8) {
256                      dj[iColumn] *= 1.0e6;
257                    } else {
258                      bool fix=false;
259                      for (j = columnStart[iColumn];
260                           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
261                        int iRow = row[j];
262                        if (bestDj[iRow]<1.0e25) {
263                          fix=true;
264                          break;
265                        }
266                      }
267                      if (fix) {
268                        dj[iColumn] *= 1.0e6;
269                      }
270                    }
271                  } else {
272                    dj[iColumn] *= 0.000001;
273                  }
[1564]274                }
275              }
276            }
277          }
278          if (fixSets)
279            djTolerance = 1.0e30;
280        }
281        delete basis;
282        delete [] sort;
283        delete [] bestDj;
284      }
[1499]285    }
[1569]286    // Do dj to get right number
287    if (type==4||type==6||type>7) {
288      double * sort = new double [numberColumns];
289      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
290        sort[iColumn]=1.0e30;
291        if (colUpper[iColumn]>colLower[iColumn]) {
292          sort[iColumn] = fabs(dj[iColumn]);
293        }
294      }
295      std::sort(sort,sort+numberColumns);
296      int last = static_cast<int>(numberColumns*fractionSmall_);
297      djTolerance = CoinMax(sort[last],1.0e-5);
298      delete [] sort;
299    }
[1564]300   
[1368]301    double primalTolerance;
302    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
303
304    int i;
305    int numberTightened = 0;
306    int numberAtBound = 0;
307    int numberContinuous = numberColumns - numberIntegers;
308
309    for (i = 0; i < numberIntegers; i++) {
310        int iColumn = integerVariable[i];
311        double value = currentSolution[iColumn];
312        double lower = colLower[iColumn];
313        double upper = colUpper[iColumn];
314        value = CoinMax(value, lower);
315        value = CoinMin(value, upper);
[1499]316        double djValue=dj[iColumn]*direction;
[1368]317#define RENS_FIX_ONLY_LOWER
318#ifndef RENS_FIX_ONLY_LOWER
319        if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
320            value = floor(value + 0.5);
321            if (value == lower || value == upper)
322                numberAtBound++;
323            newSolver->setColLower(iColumn, value);
324            newSolver->setColUpper(iColumn, value);
325            numberFixed++;
326        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
327            numberTightened++;
328            newSolver->setColLower(iColumn, floor(value));
329            newSolver->setColUpper(iColumn, ceil(value));
330        }
331#else
332        if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
[1499]333                floor(value + 0.5) == lower &&
334            djValue > djTolerance ) {
335          value = floor(value + 0.5);
336          numberAtBound++;
337          newSolver->setColLower(iColumn, value);
338          newSolver->setColUpper(iColumn, value);
339          numberFixed++;
340        } else if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
341                floor(value + 0.5) == upper &&
342                   -djValue > djTolerance && (djTolerance > 0.0||type==2)) {
343          value = floor(value + 0.5);
344          numberAtBound++;
345          newSolver->setColLower(iColumn, value);
346          newSolver->setColUpper(iColumn, value);
347          numberFixed++;
348        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0 &&
349                   djTolerance <0.0) {
[1368]350            numberTightened++;
351            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
352                value = floor(value + 0.5);
353                if (value < upper) {
354                    newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
355                    newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
356                } else {
357                    newSolver->setColLower(iColumn, upper - 1.0);
358                }
359            } else {
360                newSolver->setColLower(iColumn, floor(value));
361                newSolver->setColUpper(iColumn, ceil(value));
362            }
363        }
364#endif
365    }
[1564]366    delete [] dj;
[1368]367    if (numberFixed > numberIntegers / 5) {
[1499]368        if ( numberFixed < numberColumns / 5) {
[1368]369#define RENS_FIX_CONTINUOUS
370#ifdef RENS_FIX_CONTINUOUS
371            const double * colLower = newSolver->getColLower();
372            //const double * colUpper = newSolver->getColUpper();
373            int nAtLb = 0;
374            double sumDj = 0.0;
375            const double * dj = newSolver->getReducedCost();
376            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
377                if (!newSolver->isInteger(iColumn)) {
378                    double value = currentSolution[iColumn];
379                    if (value < colLower[iColumn] + 1.0e-8) {
380                        double djValue = dj[iColumn] * direction;
381                        nAtLb++;
382                        sumDj += djValue;
383                    }
384                }
385            }
386            if (nAtLb) {
387                // fix some continuous
388                double * sort = new double[nAtLb];
389                int * which = new int [nAtLb];
390                double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
391                int nFix2 = 0;
392                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
393                    if (!newSolver->isInteger(iColumn)) {
394                        double value = currentSolution[iColumn];
395                        if (value < colLower[iColumn] + 1.0e-8) {
396                            double djValue = dj[iColumn] * direction;
397                            if (djValue > threshold) {
398                                sort[nFix2] = -djValue;
399                                which[nFix2++] = iColumn;
400                            }
401                        }
402                    }
403                }
404                CoinSort_2(sort, sort + nFix2, which);
405                nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
406                for (int i = 0; i < nFix2; i++) {
407                    int iColumn = which[i];
408                    newSolver->setColUpper(iColumn, colLower[iColumn]);
409                }
410                delete [] sort;
411                delete [] which;
412#ifdef CLP_INVESTIGATE2
413                printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
414                       numberFixed, numberTightened, numberAtBound, nFix2);
415#endif
416            }
417#endif
418        }
419#ifdef COIN_DEVELOP
420        printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
421#endif
422        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
423                                         model_->getCutoff(), "CbcHeuristicRENS");
424        if (returnCode < 0) {
425            returnCode = 0; // returned on size
426#ifdef RENS_FIX_CONTINUOUS
427            if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
428                const double * colLower = newSolver->getColLower();
429                //const double * colUpper = newSolver->getColUpper();
430                int nAtLb = 0;
431                double sumDj = 0.0;
432                const double * dj = newSolver->getReducedCost();
433                double direction = newSolver->getObjSense();
434                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
435                    if (!newSolver->isInteger(iColumn)) {
436                        double value = currentSolution[iColumn];
437                        if (value < colLower[iColumn] + 1.0e-8) {
438                            double djValue = dj[iColumn] * direction;
439                            nAtLb++;
440                            sumDj += djValue;
441                        }
442                    }
443                }
444                if (nAtLb) {
445                    // fix some continuous
446                    double * sort = new double[nAtLb];
447                    int * which = new int [nAtLb];
448                    double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
449                    int nFix2 = 0;
450                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
451                        if (!newSolver->isInteger(iColumn)) {
452                            double value = currentSolution[iColumn];
453                            if (value < colLower[iColumn] + 1.0e-8) {
454                                double djValue = dj[iColumn] * direction;
455                                if (djValue > threshold) {
456                                    sort[nFix2] = -djValue;
457                                    which[nFix2++] = iColumn;
458                                }
459                            }
460                        }
461                    }
462                    CoinSort_2(sort, sort + nFix2, which);
463                    nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
464                    for (int i = 0; i < nFix2; i++) {
465                        int iColumn = which[i];
466                        newSolver->setColUpper(iColumn, colLower[iColumn]);
467                    }
468                    delete [] sort;
469                    delete [] which;
470#ifdef CLP_INVESTIGATE2
471                    printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
472                           numberFixed, numberTightened, numberAtBound, nFix2);
473#endif
474                }
475                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
476                                                 model_->getCutoff(), "CbcHeuristicRENS");
477#endif
478            }
479        }
480        //printf("return code %d",returnCode);
481        if ((returnCode&2) != 0) {
482            // could add cut
483            returnCode &= ~2;
484#ifdef COIN_DEVELOP
485            if (!numberTightened && numberFixed == numberAtBound)
486                printf("could add cut with %d elements\n", numberFixed);
487#endif
488        } else {
489            //printf("\n");
490        }
491    }
492
493    delete newSolver;
494    return returnCode;
495}
496// update model
497void CbcHeuristicRENS::setModel(CbcModel * model)
498{
499    model_ = model;
500}
[1432]501
Note: See TracBrowser for help on using the repository browser.