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

Last change on this file since 1564 was 1564, checked in by forrest, 9 years ago

add some heuristic variants

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