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

Last change on this file since 1569 was 1569, checked in by forrest, 8 years ago

more sos heuristics

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