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

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

add some sos heuristics

File size: 24.0 KB
Line 
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
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"
21#include "CoinWarmStartBasis.hpp"
22#include "CoinSort.hpp"
23#include "CbcBranchActual.hpp"
24#include "CbcStrategy.hpp"
25#include "CglPreProcess.hpp"
26
27// Default Constructor
28CbcHeuristicRENS::CbcHeuristicRENS()
29        : CbcHeuristic()
30{
31    numberTries_ = 0;
32    rensType_ = 0;
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;
42    rensType_ = 0;
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_;
65        rensType_ = rhs.rensType_;
66    }
67    return *this;
68}
69
70// Copy constructor
71CbcHeuristicRENS::CbcHeuristicRENS(const CbcHeuristicRENS & rhs)
72        :
73        CbcHeuristic(rhs),
74        numberTries_(rhs.numberTries_),
75        rensType_(rhs.rensType_)
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();
89    if ((numberTries_&&(rensType_&16)==0) || numberTries_>1 || (when() < 2 && bestSolution))
90        return 0;
91    numberTries_++;
92    double saveFractionSmall=fractionSmall_;
93    OsiSolverInterface * solver = model_->solver();
94
95    int numberIntegers = model_->numberIntegers();
96    const int * integerVariable = model_->integerVariable();
97
98    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
99    const double * currentSolution = newSolver->getColSolution();
100    int type = rensType_&15;
101    if (type<11)
102      newSolver->resolve();
103    double direction = newSolver->getObjSense();
104    double cutoff=model_->getCutoff();
105    newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
106    //cutoff *= direction;
107    double gap = cutoff - newSolver->getObjValue() * direction ;
108    double tolerance;
109    newSolver->getDblParam(OsiDualTolerance, tolerance) ;
110    if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<11) {
111      gap += 100.0 * tolerance;
112      int nFix = newSolver->reducedCostFix(gap);
113      if (nFix) {
114        char line [200];
115        sprintf(line, "Reduced cost fixing fixed %d variables", nFix);
116        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
117          << line
118          << CoinMessageEol;
119      }
120    } else if (type<11) {
121      return 0; // finished?
122    }
123    int numberColumns = solver->getNumCols();
124    double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
125    double djTolerance = (type!=1) ? -1.0e30 : 1.0e-4;
126    const double * colLower = newSolver->getColLower();
127    const double * colUpper = newSolver->getColUpper();
128    double * contribution = NULL;
129    int numberFixed = 0;
130    if (type==3) {
131      double total=0.0;
132      int n=0;
133      CoinWarmStartBasis * basis =
134        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
135      if (basis&&basis->getNumArtificial()) {
136        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
137          if (colUpper[iColumn]>colLower[iColumn]&&
138              basis->getStructStatus(iColumn) !=
139              CoinWarmStartBasis::basic) {
140            n++;
141            total += fabs(dj[iColumn]);
142          }
143        }
144        if (n)
145          djTolerance = (0.01*total)/static_cast<double>(n);
146        delete basis;
147      }
148    } else if (type>=5&&type<=11) {
149      /* 5 fix sets at one
150         6 fix on dj but leave unfixed SOS slacks
151         7 fix sets at one but use pi
152         8 fix all at zero but leave unfixed SOS slacks
153         9 as 8 but only fix all at zero if just one in set nonzero
154         10 fix all "stable" ones
155         11 layered approach
156      */
157      // SOS type fixing
158      bool fixSets = (type==5)||(type==7)||(type==10);
159      CoinWarmStartBasis * basis =
160        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
161      if (basis&&basis->getNumArtificial()) {
162        //const double * rowLower = solver->getRowLower();
163        const double * rowUpper = solver->getRowUpper();
164       
165        int numberRows = solver->getNumRows();
166        // Column copy
167        const CoinPackedMatrix * matrix = solver->getMatrixByCol();
168        const double * element = matrix->getElements();
169        const int * row = matrix->getIndices();
170        const CoinBigIndex * columnStart = matrix->getVectorStarts();
171        const int * columnLength = matrix->getVectorLengths();
172        double * bestDj = new double [numberRows];
173        for (int i=0;i<numberRows;i++) {
174          if (rowUpper[i]==1.0)
175            bestDj[i]=1.0e20;
176          else
177            bestDj[i]=1.0e30;
178        }
179        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
180          if (colUpper[iColumn]>colLower[iColumn]) {
181            CoinBigIndex j;
182            if (currentSolution[iColumn]>1.0e-6&&
183                currentSolution[iColumn]<0.999999) {
184              for (j = columnStart[iColumn];
185                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
186                int iRow = row[j];
187                if (bestDj[iRow]<1.0e30) {
188                  if (element[j] != 1.0)
189                    bestDj[iRow]=1.0e30;
190                  else
191                    bestDj[iRow]=1.0e25;
192                }
193              }
194            } else if ( basis->getStructStatus(iColumn) !=
195              CoinWarmStartBasis::basic) {
196              for (j = columnStart[iColumn];
197                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
198                int iRow = row[j];
199                if (bestDj[iRow]<1.0e25) {
200                  if (element[j] != 1.0)
201                    bestDj[iRow]=1.0e30;
202                  else
203                    bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]);
204                }
205              }
206            }
207          }
208        }
209        // Just leave one slack in each set
210        {
211          const double * objective = newSolver->getObjCoefficients();
212          int * best = new int [numberRows];
213          double * cheapest = new double[numberRows];
214          for (int i=0;i<numberRows;i++) {
215            best[i]=-1;
216            cheapest[i]=COIN_DBL_MAX;
217          }
218          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
219            if (colUpper[iColumn]>colLower[iColumn]) {
220              if (columnLength[iColumn]==1) {
221                CoinBigIndex j = columnStart[iColumn];
222                int iRow = row[j];
223                if (bestDj[iRow]<1.0e30) {
224                  double obj = direction*objective[iColumn];
225                  if (obj<cheapest[iRow]) {
226                    cheapest[iRow]=obj;
227                    best[iRow]=iColumn;
228                  }
229                }
230              }
231            }
232          }
233          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
234            if (colUpper[iColumn]>colLower[iColumn]) {
235              if (columnLength[iColumn]==1) {
236                CoinBigIndex j = columnStart[iColumn];
237                int iRow = row[j];
238                if (bestDj[iRow]<1.0e30) {
239                  if (best[iRow]!=-1&&iColumn!=best[iRow]) {
240                    newSolver->setColUpper(iColumn,0.0);
241                  }
242                }
243              }
244            }
245          }
246          delete [] best;
247          delete [] cheapest;
248        }
249        int nSOS=0;
250        double * sort = new double [numberRows];
251        const double * pi = newSolver->getRowPrice();
252        if (type==11) {
253          contribution = new double [numberRows];
254          for (int i=0;i<numberRows;i++) {
255            if (bestDj[i]<1.0e30) 
256              contribution[i]=0.0;
257            else
258              contribution[i]=-1.0;
259          }
260        }
261        for (int i=0;i<numberRows;i++) {
262          if (bestDj[i]<1.0e30) {
263            if (type==5)
264              sort[nSOS++]=bestDj[i];
265            else if (type==7)
266              sort[nSOS++]=-fabs(pi[i]);
267            else
268              sort[nSOS++]=fabs(pi[i]);
269          }
270        }
271        if (10*nSOS>8*numberRows) {
272          if (type<10) {
273            std::sort(sort,sort+nSOS);
274            int last = static_cast<int>(nSOS*0.9*fractionSmall_);
275            double tolerance = sort[last];
276            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
277              if (colUpper[iColumn]>colLower[iColumn]) {
278                CoinBigIndex j;
279                if (currentSolution[iColumn]<=1.0e-6||
280                    currentSolution[iColumn]>=0.999999) {
281                  if (fixSets) {
282                    for (j = columnStart[iColumn];
283                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
284                      int iRow = row[j];
285                      double useDj;
286                      if (type==5) 
287                        useDj = bestDj[iRow];
288                      else if (type==7)
289                        useDj= -fabs(pi[iRow]);
290                      else
291                        useDj= fabs(pi[iRow]);
292                      if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
293                        numberFixed++;
294                        if (currentSolution[iColumn]<=1.0e-6)
295                          newSolver->setColUpper(iColumn,0.0);
296                        else if (currentSolution[iColumn]>=0.999999) 
297                          newSolver->setColLower(iColumn,1.0);
298                      }
299                    }
300                  } else if (columnLength[iColumn]==1) {
301                    // leave more slacks
302                    int iRow = row[columnStart[iColumn]];
303                    if (bestDj[iRow]<1.0e30) {
304                      // fake dj
305                      dj[iColumn] *= 0.000001;
306                    }
307                  } else if (type==8||type==9) {
308                    if (currentSolution[iColumn]<=1.0e-6) {
309                      if (type==8) {
310                        dj[iColumn] *= 1.0e6;
311                      } else {
312                        bool fix=false;
313                        for (j = columnStart[iColumn];
314                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
315                          int iRow = row[j];
316                          if (bestDj[iRow]<1.0e25) {
317                            fix=true;
318                            break;
319                          }
320                        }
321                        if (fix) {
322                          dj[iColumn] *= 1.0e6;
323                        }
324                      }
325                    } else {
326                      dj[iColumn] *= 0.000001;
327                    }
328                  }
329                }
330              }
331            }
332            if (fixSets)
333              djTolerance = 1.0e30;
334          } else if (type==10) {
335            double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
336            char * mark = new char [numberColumns];
337            char * nonzero = new char [numberColumns];
338            double factor=CoinMax(1.000001,fractionSmall_);
339            fractionSmall_ = 0.5;
340            // loosen up
341            for (int i=0;i<numberRows;i++) {
342              if (bestDj[i]>=1.0e30) {
343                newSolver->setRowUpper(i,factor*saveUpper[i]);
344              }
345            }
346            newSolver->resolve();
347            const double * solution = newSolver->getColSolution();
348            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
349              mark[iColumn]=0;
350              nonzero[iColumn]=0;
351              if (colUpper[iColumn]>colLower[iColumn]&&
352                  solution[iColumn]>0.9999)
353                mark[iColumn]=1;
354              else if (solution[iColumn]>0.00001)
355                nonzero[iColumn]=1;
356            }
357            // slightly small
358            for (int i=0;i<numberRows;i++) {
359              if (bestDj[i]>=1.0e30) {
360                newSolver->setRowUpper(i,saveUpper[i]*0.9999);
361              }
362            }
363            newSolver->resolve();
364            int nCheck=2;
365            if (newSolver->isProvenOptimal()) {
366              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
367                if (colUpper[iColumn]>colLower[iColumn]&&
368                    solution[iColumn]>0.9999)
369                  mark[iColumn]++;
370                else if (solution[iColumn]>0.00001)
371                  nonzero[iColumn]=1;
372              }
373            } else {
374              nCheck=1;
375            }
376            // correct values
377            for (int i=0;i<numberRows;i++) {
378              if (bestDj[i]>=1.0e30) {
379                newSolver->setRowUpper(i,saveUpper[i]);
380              }
381            }
382            newSolver->resolve();
383            int nFixed=0;
384            int nFixedToZero=0;
385            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
386              if (colUpper[iColumn]>colLower[iColumn]) {
387                if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) {
388                  newSolver->setColLower(iColumn,1.0);
389                  nFixed++;
390                } else if (!mark[iColumn]&&!nonzero[iColumn]&&
391                           columnLength[iColumn]>1&&solution[iColumn]<0.00001) {
392                  newSolver->setColUpper(iColumn,0.0);
393                  nFixedToZero++;
394                }
395              }
396            }
397            char line[100];
398            sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)",
399                    heuristicName(),
400                    nFixed,nFixedToZero);
401            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
402              << line
403              << CoinMessageEol;
404            delete [] mark;
405            delete []nonzero;
406            delete [] saveUpper;
407            numberFixed=numberColumns;
408            djTolerance = 1.0e30;
409          }
410        }
411        delete basis;
412        delete [] sort;
413        delete [] bestDj;
414        if (10*nSOS<=8*numberRows) {
415          // give up
416          delete [] contribution;
417          delete newSolver;
418          return 0;
419        }
420      }
421    }
422    // Do dj to get right number
423    if (type==4||type==6||(type>7&&type<10)) {
424      double * sort = new double [numberColumns];
425      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
426        sort[iColumn]=1.0e30;
427        if (colUpper[iColumn]>colLower[iColumn]) {
428          sort[iColumn] = fabs(dj[iColumn]);
429        }
430      }
431      std::sort(sort,sort+numberColumns);
432      int last = static_cast<int>(numberColumns*fractionSmall_);
433      djTolerance = CoinMax(sort[last],1.0e-5);
434      delete [] sort;
435    } else if (type==11) {
436      // Do layered in a different way
437      int numberRows = solver->getNumRows();
438      // Column copy
439      const CoinPackedMatrix * matrix = newSolver->getMatrixByCol();
440      const double * element = matrix->getElements();
441      const int * row = matrix->getIndices();
442      const CoinBigIndex * columnStart = matrix->getVectorStarts();
443      const int * columnLength = matrix->getVectorLengths();
444      int * whichRow = new int[numberRows];
445      int * whichSet = new int [numberColumns];
446      int nSOS=0;
447      for (int i=0;i<numberRows;i++) {
448        whichRow[i]=0;
449        if (!contribution[i])
450          nSOS++;
451      }
452      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
453        whichSet[iColumn]=-2;
454        if (colUpper[iColumn]>colLower[iColumn]) {
455          CoinBigIndex j;
456          double sum=0.0;
457          int iSOS=-1;
458          int n=0;
459          for (j = columnStart[iColumn];
460               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
461            int iRow = row[j];
462            if (contribution[iRow]>=0.0) {
463              iSOS=iRow;
464              n++;
465            } else {
466              sum += fabs(element[j]);
467            }
468          }
469          if (n>1)
470            printf("Too many SOS entries (%d) for column %d\n",
471                   n,iColumn);
472          if (sum) {
473            assert (iSOS>=0);
474            contribution[iSOS] += sum;
475            whichRow[iSOS]++;
476            whichSet[iColumn]=iSOS;
477          } else {
478            whichSet[iColumn]=iSOS+numberRows;
479          }
480        }
481      }
482      int * chunk = new int [numberRows];
483      for (int i=0;i<numberRows;i++) {
484        chunk[i]=-1;
485        if (whichRow[i]) {
486          contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]);
487        } else {
488          contribution[i] = COIN_DBL_MAX;
489        }
490        whichRow[i]=i;
491      }
492      newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
493      double * saveLower = CoinCopyOfArray(colLower,numberColumns);
494      double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
495      CoinSort_2(contribution,contribution+numberRows,whichRow);
496      // Set do nothing solution
497      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
498        if(whichSet[iColumn]>=numberRows)
499          newSolver->setColLower(iColumn,1.0);
500      }
501      newSolver->resolve();
502      int nChunk = (nSOS+9)/10;
503      int nPass=0;
504      int inChunk=0;
505      for (int i=0;i<nSOS;i++) {
506        chunk[whichRow[i]]=nPass;
507        inChunk++;
508        if (inChunk==nChunk) {
509          inChunk=0;
510          // last two together
511          if (i+nChunk<nSOS)
512            nPass++;
513        }
514      }
515      // adjust
516      nPass++;
517      for (int iPass=0;iPass<nPass;iPass++) {
518        // fix last chunk and unfix this chunk
519        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
520          int iSOS = whichSet[iColumn];
521          if (iSOS>=0) {
522            if (iSOS>=numberRows)
523              iSOS-=numberRows;
524            if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) {
525              newSolver->setColLower(iColumn,1.0);
526            } else if (chunk[iSOS]==iPass) {
527              newSolver->setColLower(iColumn,saveLower[iColumn]);
528              newSolver->setColUpper(iColumn,saveUpper[iColumn]);
529            }
530          }
531        }
532        // solve
533        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
534                                         model_->getCutoff(), "CbcHeuristicRENS");
535        if (returnCode < 0) {
536            returnCode = 0; // returned on size
537            break;
538        } else if ((returnCode&1)==0) {
539          // no good
540          break;
541        }
542      }
543      if ((returnCode&2) != 0) {
544        // could add cut
545        returnCode &= ~2;
546      }
547      delete [] chunk;
548      delete [] saveLower;
549      delete [] saveUpper;
550      delete [] whichRow;
551      delete [] whichSet;
552      delete [] contribution;
553      delete newSolver;
554      return returnCode;
555    }
556   
557    double primalTolerance;
558    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
559
560    int i;
561    int numberTightened = 0;
562    int numberAtBound = 0;
563    int numberContinuous = numberColumns - numberIntegers;
564
565    for (i = 0; i < numberIntegers; i++) {
566        int iColumn = integerVariable[i];
567        double value = currentSolution[iColumn];
568        double lower = colLower[iColumn];
569        double upper = colUpper[iColumn];
570        value = CoinMax(value, lower);
571        value = CoinMin(value, upper);
572        double djValue=dj[iColumn]*direction;
573#define RENS_FIX_ONLY_LOWER
574#ifndef RENS_FIX_ONLY_LOWER
575        if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
576            value = floor(value + 0.5);
577            if (value == lower || value == upper)
578                numberAtBound++;
579            newSolver->setColLower(iColumn, value);
580            newSolver->setColUpper(iColumn, value);
581            numberFixed++;
582        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
583            numberTightened++;
584            newSolver->setColLower(iColumn, floor(value));
585            newSolver->setColUpper(iColumn, ceil(value));
586        }
587#else
588        if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
589                floor(value + 0.5) == lower &&
590            djValue > djTolerance ) {
591          value = floor(value + 0.5);
592          numberAtBound++;
593          newSolver->setColLower(iColumn, value);
594          newSolver->setColUpper(iColumn, value);
595          numberFixed++;
596        } else if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
597                floor(value + 0.5) == upper &&
598                   -djValue > djTolerance && (djTolerance > 0.0||type==2)) {
599          value = floor(value + 0.5);
600          numberAtBound++;
601          newSolver->setColLower(iColumn, value);
602          newSolver->setColUpper(iColumn, value);
603          numberFixed++;
604        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0 &&
605                   djTolerance <0.0) {
606            numberTightened++;
607            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
608                value = floor(value + 0.5);
609                if (value < upper) {
610                    newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
611                    newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
612                } else {
613                    newSolver->setColLower(iColumn, upper - 1.0);
614                }
615            } else {
616                newSolver->setColLower(iColumn, floor(value));
617                newSolver->setColUpper(iColumn, ceil(value));
618            }
619        }
620#endif
621    }
622    delete [] dj;
623    if (numberFixed > numberIntegers / 5) {
624        if ( numberFixed < numberColumns / 5) {
625#define RENS_FIX_CONTINUOUS
626#ifdef RENS_FIX_CONTINUOUS
627            const double * colLower = newSolver->getColLower();
628            //const double * colUpper = newSolver->getColUpper();
629            int nAtLb = 0;
630            double sumDj = 0.0;
631            const double * dj = newSolver->getReducedCost();
632            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
633                if (!newSolver->isInteger(iColumn)) {
634                    double value = currentSolution[iColumn];
635                    if (value < colLower[iColumn] + 1.0e-8) {
636                        double djValue = dj[iColumn] * direction;
637                        nAtLb++;
638                        sumDj += djValue;
639                    }
640                }
641            }
642            if (nAtLb) {
643                // fix some continuous
644                double * sort = new double[nAtLb];
645                int * which = new int [nAtLb];
646                double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
647                int nFix2 = 0;
648                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
649                    if (!newSolver->isInteger(iColumn)) {
650                        double value = currentSolution[iColumn];
651                        if (value < colLower[iColumn] + 1.0e-8) {
652                            double djValue = dj[iColumn] * direction;
653                            if (djValue > threshold) {
654                                sort[nFix2] = -djValue;
655                                which[nFix2++] = iColumn;
656                            }
657                        }
658                    }
659                }
660                CoinSort_2(sort, sort + nFix2, which);
661                nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
662                for (int i = 0; i < nFix2; i++) {
663                    int iColumn = which[i];
664                    newSolver->setColUpper(iColumn, colLower[iColumn]);
665                }
666                delete [] sort;
667                delete [] which;
668#ifdef CLP_INVESTIGATE2
669                printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
670                       numberFixed, numberTightened, numberAtBound, nFix2);
671#endif
672            }
673#endif
674        }
675#ifdef COIN_DEVELOP
676        printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
677#endif
678        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
679                                         model_->getCutoff(), "CbcHeuristicRENS");
680        if (returnCode < 0) {
681            returnCode = 0; // returned on size
682#ifdef RENS_FIX_CONTINUOUS
683            if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
684                const double * colLower = newSolver->getColLower();
685                //const double * colUpper = newSolver->getColUpper();
686                int nAtLb = 0;
687                double sumDj = 0.0;
688                const double * dj = newSolver->getReducedCost();
689                double direction = newSolver->getObjSense();
690                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
691                    if (!newSolver->isInteger(iColumn)) {
692                        double value = currentSolution[iColumn];
693                        if (value < colLower[iColumn] + 1.0e-8) {
694                            double djValue = dj[iColumn] * direction;
695                            nAtLb++;
696                            sumDj += djValue;
697                        }
698                    }
699                }
700                if (nAtLb) {
701                    // fix some continuous
702                    double * sort = new double[nAtLb];
703                    int * which = new int [nAtLb];
704                    double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
705                    int nFix2 = 0;
706                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
707                        if (!newSolver->isInteger(iColumn)) {
708                            double value = currentSolution[iColumn];
709                            if (value < colLower[iColumn] + 1.0e-8) {
710                                double djValue = dj[iColumn] * direction;
711                                if (djValue > threshold) {
712                                    sort[nFix2] = -djValue;
713                                    which[nFix2++] = iColumn;
714                                }
715                            }
716                        }
717                    }
718                    CoinSort_2(sort, sort + nFix2, which);
719                    nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
720                    for (int i = 0; i < nFix2; i++) {
721                        int iColumn = which[i];
722                        newSolver->setColUpper(iColumn, colLower[iColumn]);
723                    }
724                    delete [] sort;
725                    delete [] which;
726#ifdef CLP_INVESTIGATE2
727                    printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
728                           numberFixed, numberTightened, numberAtBound, nFix2);
729#endif
730                }
731                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
732                                                 model_->getCutoff(), "CbcHeuristicRENS");
733#endif
734            }
735        }
736        //printf("return code %d",returnCode);
737        if ((returnCode&2) != 0) {
738            // could add cut
739            returnCode &= ~2;
740#ifdef COIN_DEVELOP
741            if (!numberTightened && numberFixed == numberAtBound)
742                printf("could add cut with %d elements\n", numberFixed);
743#endif
744        } else {
745            //printf("\n");
746        }
747    }
748    //delete [] whichRow;
749    //delete [] contribution;
750    delete newSolver;
751    fractionSmall_ = saveFractionSmall;
752    return returnCode;
753}
754// update model
755void CbcHeuristicRENS::setModel(CbcModel * model)
756{
757    model_ = model;
758}
759
Note: See TracBrowser for help on using the repository browser.