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

Last change on this file since 1591 was 1591, checked in by forrest, 10 years ago

modifications to heuristics and allow missing out some printout

File size: 29.4 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<12)
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<12) {
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<12) {
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<=12) {
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 fix all "stable" ones - approach 2
156         12 layered approach
157      */
158      // SOS type fixing
159      bool fixSets = (type==5)||(type==7)||(type==10)||(type==11);
160      CoinWarmStartBasis * basis =
161        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
162      if (basis&&basis->getNumArtificial()) {
163        //const double * rowLower = solver->getRowLower();
164        const double * rowUpper = solver->getRowUpper();
165       
166        int numberRows = solver->getNumRows();
167        // Column copy
168        const CoinPackedMatrix * matrix = solver->getMatrixByCol();
169        const double * element = matrix->getElements();
170        const int * row = matrix->getIndices();
171        const CoinBigIndex * columnStart = matrix->getVectorStarts();
172        const int * columnLength = matrix->getVectorLengths();
173        double * bestDj = new double [numberRows];
174        for (int i=0;i<numberRows;i++) {
175          if (rowUpper[i]==1.0)
176            bestDj[i]=1.0e20;
177          else
178            bestDj[i]=1.0e30;
179        }
180        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
181          if (colUpper[iColumn]>colLower[iColumn]) {
182            CoinBigIndex j;
183            if (currentSolution[iColumn]>1.0e-6&&
184                currentSolution[iColumn]<0.999999) {
185              for (j = columnStart[iColumn];
186                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
187                int iRow = row[j];
188                if (bestDj[iRow]<1.0e30) {
189                  if (element[j] != 1.0)
190                    bestDj[iRow]=1.0e30;
191                  else
192                    bestDj[iRow]=1.0e25;
193                }
194              }
195            } else if ( basis->getStructStatus(iColumn) !=
196              CoinWarmStartBasis::basic) {
197              for (j = columnStart[iColumn];
198                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
199                int iRow = row[j];
200                if (bestDj[iRow]<1.0e25) {
201                  if (element[j] != 1.0)
202                    bestDj[iRow]=1.0e30;
203                  else
204                    bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]);
205                }
206              }
207            }
208          }
209        }
210        // Just leave one slack in each set
211        {
212          const double * objective = newSolver->getObjCoefficients();
213          int * best = new int [numberRows];
214          double * cheapest = new double[numberRows];
215          for (int i=0;i<numberRows;i++) {
216            best[i]=-1;
217            cheapest[i]=COIN_DBL_MAX;
218          }
219          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
220            if (colUpper[iColumn]>colLower[iColumn]) {
221              if (columnLength[iColumn]==1) {
222                CoinBigIndex j = columnStart[iColumn];
223                int iRow = row[j];
224                if (bestDj[iRow]<1.0e30) {
225                  double obj = direction*objective[iColumn];
226                  if (obj<cheapest[iRow]) {
227                    cheapest[iRow]=obj;
228                    best[iRow]=iColumn;
229                  }
230                }
231              }
232            }
233          }
234          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
235            if (colUpper[iColumn]>colLower[iColumn]) {
236              if (columnLength[iColumn]==1) {
237                CoinBigIndex j = columnStart[iColumn];
238                int iRow = row[j];
239                if (bestDj[iRow]<1.0e30) {
240                  if (best[iRow]!=-1&&iColumn!=best[iRow]) {
241                    newSolver->setColUpper(iColumn,0.0);
242                  }
243                }
244              }
245            }
246          }
247          delete [] best;
248          delete [] cheapest;
249        }
250        int nSOS=0;
251        double * sort = new double [numberRows];
252        const double * pi = newSolver->getRowPrice();
253        if (type==12) {
254          contribution = new double [numberRows];
255          for (int i=0;i<numberRows;i++) {
256            if (bestDj[i]<1.0e30) 
257              contribution[i]=0.0;
258            else
259              contribution[i]=-1.0;
260          }
261        }
262        for (int i=0;i<numberRows;i++) {
263          if (bestDj[i]<1.0e30) {
264            if (type==5)
265              sort[nSOS++]=bestDj[i];
266            else if (type==7)
267              sort[nSOS++]=-fabs(pi[i]);
268            else
269              sort[nSOS++]=fabs(pi[i]);
270          }
271        }
272        if (10*nSOS>8*numberRows) {
273          if (type<10) {
274            std::sort(sort,sort+nSOS);
275            int last = static_cast<int>(nSOS*0.9*fractionSmall_);
276            double tolerance = sort[last];
277            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
278              if (colUpper[iColumn]>colLower[iColumn]) {
279                CoinBigIndex j;
280                if (currentSolution[iColumn]<=1.0e-6||
281                    currentSolution[iColumn]>=0.999999) {
282                  if (fixSets) {
283                    for (j = columnStart[iColumn];
284                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
285                      int iRow = row[j];
286                      double useDj;
287                      if (type==5) 
288                        useDj = bestDj[iRow];
289                      else if (type==7)
290                        useDj= -fabs(pi[iRow]);
291                      else
292                        useDj= fabs(pi[iRow]);
293                      if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
294                        numberFixed++;
295                        if (currentSolution[iColumn]<=1.0e-6)
296                          newSolver->setColUpper(iColumn,0.0);
297                        else if (currentSolution[iColumn]>=0.999999) 
298                          newSolver->setColLower(iColumn,1.0);
299                      }
300                    }
301                  } else if (columnLength[iColumn]==1) {
302                    // leave more slacks
303                    int iRow = row[columnStart[iColumn]];
304                    if (bestDj[iRow]<1.0e30) {
305                      // fake dj
306                      dj[iColumn] *= 0.000001;
307                    }
308                  } else if (type==8||type==9) {
309                    if (currentSolution[iColumn]<=1.0e-6) {
310                      if (type==8) {
311                        dj[iColumn] *= 1.0e6;
312                      } else {
313                        bool fix=false;
314                        for (j = columnStart[iColumn];
315                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
316                          int iRow = row[j];
317                          if (bestDj[iRow]<1.0e25) {
318                            fix=true;
319                            break;
320                          }
321                        }
322                        if (fix) {
323                          dj[iColumn] *= 1.0e6;
324                        }
325                      }
326                    } else {
327                      dj[iColumn] *= 0.000001;
328                    }
329                  }
330                }
331              }
332            }
333            if (fixSets)
334              djTolerance = 1.0e30;
335          } else if (type==10) {
336            double * saveUpper = new double [numberRows];
337            memset(saveUpper,0,numberRows*sizeof(double));
338            char * mark = new char [numberColumns];
339            char * nonzero = new char [numberColumns];
340            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
341              if (colUpper[iColumn]>colLower[iColumn]) {
342                CoinBigIndex j;
343                for (j = columnStart[iColumn];
344                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
345                  int iRow = row[j];
346                  saveUpper[iRow] += element[j];
347                }
348              }
349            }
350            double sum=0.0;
351            double sumRhs=0.0;
352            const double * rowUpper = newSolver->getRowUpper();
353            for (int i=0;i<numberRows;i++) {
354              if (bestDj[i]>=1.0e30) {
355                sum += saveUpper[i];
356                sumRhs += rowUpper[i];
357              }
358            }
359            double averagePerSet = sum/static_cast<double>(numberRows);
360            // allow this extra
361            double factor = averagePerSet*fractionSmall_*numberRows;
362            factor = 1.0+factor/sumRhs;
363            fractionSmall_ = 0.5;
364            memcpy(saveUpper,rowUpper,numberRows*sizeof(double));
365            // loosen up
366            for (int i=0;i<numberRows;i++) {
367              if (bestDj[i]>=1.0e30) {
368                newSolver->setRowUpper(i,factor*saveUpper[i]);
369              }
370            }
371            newSolver->resolve();
372            const double * solution = newSolver->getColSolution();
373            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
374              mark[iColumn]=0;
375              nonzero[iColumn]=0;
376              if (colUpper[iColumn]>colLower[iColumn]&&
377                  solution[iColumn]>0.9999)
378                mark[iColumn]=1;
379              else if (solution[iColumn]>0.00001)
380                nonzero[iColumn]=1;
381            }
382            // slightly small
383            for (int i=0;i<numberRows;i++) {
384              if (bestDj[i]>=1.0e30) {
385                newSolver->setRowUpper(i,saveUpper[i]*0.9999);
386              }
387            }
388            newSolver->resolve();
389            int nCheck=2;
390            if (newSolver->isProvenOptimal()) {
391              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
392                if (colUpper[iColumn]>colLower[iColumn]&&
393                    solution[iColumn]>0.9999)
394                  mark[iColumn]++;
395                else if (solution[iColumn]>0.00001)
396                  nonzero[iColumn]=1;
397              }
398            } else {
399              nCheck=1;
400            }
401            // correct values
402            for (int i=0;i<numberRows;i++) {
403              if (bestDj[i]>=1.0e30) {
404                newSolver->setRowUpper(i,saveUpper[i]);
405              }
406            }
407            newSolver->resolve();
408            int nFixed=0;
409            int nFixedToZero=0;
410            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
411              if (colUpper[iColumn]>colLower[iColumn]) {
412                if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) {
413                  newSolver->setColLower(iColumn,1.0);
414                  nFixed++;
415                } else if (!mark[iColumn]&&!nonzero[iColumn]&&
416                           columnLength[iColumn]>1&&solution[iColumn]<0.00001) {
417                  newSolver->setColUpper(iColumn,0.0);
418                  nFixedToZero++;
419                }
420              }
421            }
422            char line[100];
423            sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)",
424                    heuristicName(),
425                    nFixed,nFixedToZero);
426            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
427              << line
428              << CoinMessageEol;
429            delete [] mark;
430            delete []nonzero;
431            delete [] saveUpper;
432            numberFixed=numberColumns;
433            djTolerance = 1.0e30;
434          } else if (type==11) {
435            double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
436            char * mark = new char [numberColumns];
437            char * nonzero = new char [numberColumns];
438            // save basis and solution
439            CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(newSolver->getWarmStart()) ;
440            assert(basis != NULL);
441            double * saveSolution = 
442              CoinCopyOfArray(newSolver->getColSolution(), 
443                              numberColumns);
444            double factors[] = {1.1,1.05,1.01,0.98};
445            int nPass = (sizeof(factors)/sizeof(double))-1;
446            double factor=factors[0];
447            double proportion = fractionSmall_;
448            fractionSmall_ = 0.5;
449            // loosen up
450            for (int i=0;i<numberRows;i++) {
451              if (bestDj[i]>=1.0e30) {
452                newSolver->setRowUpper(i,factor*saveUpper[i]);
453              }
454            }
455            bool takeHint;
456            OsiHintStrength strength;
457            newSolver->getHintParam(OsiDoDualInResolve, takeHint, strength);
458            newSolver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
459            newSolver->resolve();
460            newSolver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
461            const double * solution = newSolver->getColSolution();
462            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
463              mark[iColumn]=0;
464              nonzero[iColumn]=0;
465              if (colUpper[iColumn]>colLower[iColumn]&&
466                  solution[iColumn]>0.9999)
467                mark[iColumn]=1;
468              else if (solution[iColumn]>0.00001)
469                nonzero[iColumn]=1;
470            }
471            int nCheck=2;
472            for (int iPass=0;iPass<nPass;iPass++) {
473              // smaller
474              factor = factors[iPass+1];
475              for (int i=0;i<numberRows;i++) {
476                if (bestDj[i]>=1.0e30) {
477                  newSolver->setRowUpper(i,saveUpper[i]*factor);
478                }
479              }
480              newSolver->resolve();
481              if (newSolver->isProvenOptimal()) {
482                nCheck++;
483                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
484                  if (colUpper[iColumn]>colLower[iColumn]&&
485                      solution[iColumn]>0.9999)
486                    mark[iColumn]++;
487                  else if (solution[iColumn]>0.00001)
488                    nonzero[iColumn]++;
489                }
490              }
491            }
492            // correct values
493            for (int i=0;i<numberRows;i++) {
494              if (bestDj[i]>=1.0e30) {
495                newSolver->setRowUpper(i,saveUpper[i]);
496              }
497            }
498            newSolver->setColSolution(saveSolution);
499            delete [] saveSolution;
500            newSolver->setWarmStart(basis);
501            delete basis ;
502            newSolver->setHintParam(OsiDoDualInResolve, takeHint, strength);
503            newSolver->resolve();
504            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
505              if (colUpper[iColumn]>colLower[iColumn]&&
506                  solution[iColumn]>0.9999)
507                mark[iColumn]++;
508              else if (solution[iColumn]>0.00001)
509                nonzero[iColumn]++;
510            }
511            int nFixed=0;
512            int numberSetsToFix = static_cast<int>(nSOS*(1.0-proportion));
513            int * mixed = new int[numberRows];
514            memset(mixed,0,numberRows*sizeof(int));
515            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
516              if (colUpper[iColumn]>colLower[iColumn]) {
517                int iSOS=-1;
518                for (CoinBigIndex j = columnStart[iColumn];
519                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
520                  int iRow = row[j];
521                  if (bestDj[iRow]<1.0e25) {
522                    iSOS=iRow;
523                    break;
524                  }
525                }
526                if (iSOS>=0) {
527                  int numberTimesAtOne = mark[iColumn];
528                  int numberTimesNonZero = nonzero[iColumn]+
529                    numberTimesAtOne;
530                  if (numberTimesAtOne<nCheck&&
531                      numberTimesNonZero) {
532                    mixed[iSOS]+=
533                      CoinMin(numberTimesNonZero,
534                              nCheck-numberTimesNonZero);
535                  } 
536                }
537              }
538            }
539            int nFix=0;
540            for (int i=0;i<numberRows;i++) {
541              if (bestDj[i]<1.0e25) {
542                sort[nFix] = -bestDj[i]+1.0e8*mixed[i];
543                mixed[nFix++]=i;
544              }
545            }
546            CoinSort_2(sort,sort+nFix,mixed);
547            nFix = CoinMin(nFix,numberSetsToFix);
548            memset(sort,0,sizeof(double)*numberRows);
549            for (int i=0;i<nFix;i++)
550              sort[mixed[i]]=1.0;
551            delete [] mixed;
552            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
553              if (colUpper[iColumn]>colLower[iColumn]) {
554                if (solution[iColumn]>0.9999) {
555                  int iSOS=-1;
556                  for (CoinBigIndex j = columnStart[iColumn];
557                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
558                    int iRow = row[j];
559                    if (bestDj[iRow]<1.0e25) {
560                      iSOS=iRow;
561                      break;
562                    }
563                  }
564                  if (iSOS>=0&&sort[iSOS]) {
565                    newSolver->setColLower(iColumn,1.0);
566                    nFixed++;
567                  }
568                }
569              }
570            }
571            char line[100];
572            sprintf(line,"Heuristic %s fixed %d to one (%d sets)",
573                    heuristicName(),
574                    nFixed,nSOS);
575            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
576              << line
577              << CoinMessageEol;
578            delete [] mark;
579            delete [] nonzero;
580            delete [] saveUpper;
581            numberFixed=numberColumns;
582            djTolerance = 1.0e30;
583          }
584        }
585        delete basis;
586        delete [] sort;
587        delete [] bestDj;
588        if (10*nSOS<=8*numberRows) {
589          // give up
590          delete [] contribution;
591          delete newSolver;
592          return 0;
593        }
594      }
595    }
596    // Do dj to get right number
597    if (type==4||type==6||(type>7&&type<10)) {
598      double * sort = new double [numberColumns];
599      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
600        sort[iColumn]=1.0e30;
601        if (colUpper[iColumn]>colLower[iColumn]) {
602          sort[iColumn] = fabs(dj[iColumn]);
603        }
604      }
605      std::sort(sort,sort+numberColumns);
606      int last = static_cast<int>(numberColumns*fractionSmall_);
607      djTolerance = CoinMax(sort[last],1.0e-5);
608      delete [] sort;
609    } else if (type==12) {
610      // Do layered in a different way
611      int numberRows = solver->getNumRows();
612      // Column copy
613      const CoinPackedMatrix * matrix = newSolver->getMatrixByCol();
614      const double * element = matrix->getElements();
615      const int * row = matrix->getIndices();
616      const CoinBigIndex * columnStart = matrix->getVectorStarts();
617      const int * columnLength = matrix->getVectorLengths();
618      int * whichRow = new int[numberRows];
619      int * whichSet = new int [numberColumns];
620      int nSOS=0;
621      for (int i=0;i<numberRows;i++) {
622        whichRow[i]=0;
623        if (!contribution[i])
624          nSOS++;
625      }
626      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
627        whichSet[iColumn]=-2;
628        if (colUpper[iColumn]>colLower[iColumn]) {
629          CoinBigIndex j;
630          double sum=0.0;
631          int iSOS=-1;
632          int n=0;
633          for (j = columnStart[iColumn];
634               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
635            int iRow = row[j];
636            if (contribution[iRow]>=0.0) {
637              iSOS=iRow;
638              n++;
639            } else {
640              sum += fabs(element[j]);
641            }
642          }
643          if (n>1)
644            printf("Too many SOS entries (%d) for column %d\n",
645                   n,iColumn);
646          if (sum) {
647            assert (iSOS>=0);
648            contribution[iSOS] += sum;
649            whichRow[iSOS]++;
650            whichSet[iColumn]=iSOS;
651          } else {
652            whichSet[iColumn]=iSOS+numberRows;
653          }
654        }
655      }
656      int * chunk = new int [numberRows];
657      for (int i=0;i<numberRows;i++) {
658        chunk[i]=-1;
659        if (whichRow[i]) {
660          contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]);
661        } else {
662          contribution[i] = COIN_DBL_MAX;
663        }
664        whichRow[i]=i;
665      }
666      newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
667      double * saveLower = CoinCopyOfArray(colLower,numberColumns);
668      double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
669      CoinSort_2(contribution,contribution+numberRows,whichRow);
670      // Set do nothing solution
671      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
672        if(whichSet[iColumn]>=numberRows)
673          newSolver->setColLower(iColumn,1.0);
674      }
675      newSolver->resolve();
676      int nChunk = (nSOS+9)/10;
677      int nPass=0;
678      int inChunk=0;
679      for (int i=0;i<nSOS;i++) {
680        chunk[whichRow[i]]=nPass;
681        inChunk++;
682        if (inChunk==nChunk) {
683          inChunk=0;
684          // last two together
685          if (i+nChunk<nSOS)
686            nPass++;
687        }
688      }
689      // adjust
690      nPass++;
691      for (int iPass=0;iPass<nPass;iPass++) {
692        // fix last chunk and unfix this chunk
693        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
694          int iSOS = whichSet[iColumn];
695          if (iSOS>=0) {
696            if (iSOS>=numberRows)
697              iSOS-=numberRows;
698            if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) {
699              newSolver->setColLower(iColumn,1.0);
700            } else if (chunk[iSOS]==iPass) {
701              newSolver->setColLower(iColumn,saveLower[iColumn]);
702              newSolver->setColUpper(iColumn,saveUpper[iColumn]);
703            }
704          }
705        }
706        // solve
707        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
708                                         model_->getCutoff(), "CbcHeuristicRENS");
709        if (returnCode < 0) {
710            returnCode = 0; // returned on size
711            break;
712        } else if ((returnCode&1)==0) {
713          // no good
714          break;
715        }
716      }
717      if ((returnCode&2) != 0) {
718        // could add cut
719        returnCode &= ~2;
720      }
721      delete [] chunk;
722      delete [] saveLower;
723      delete [] saveUpper;
724      delete [] whichRow;
725      delete [] whichSet;
726      delete [] contribution;
727      delete newSolver;
728      return returnCode;
729    }
730   
731    double primalTolerance;
732    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
733
734    int i;
735    int numberTightened = 0;
736    int numberAtBound = 0;
737    int numberContinuous = numberColumns - numberIntegers;
738
739    for (i = 0; i < numberIntegers; i++) {
740        int iColumn = integerVariable[i];
741        double value = currentSolution[iColumn];
742        double lower = colLower[iColumn];
743        double upper = colUpper[iColumn];
744        value = CoinMax(value, lower);
745        value = CoinMin(value, upper);
746        double djValue=dj[iColumn]*direction;
747#define RENS_FIX_ONLY_LOWER
748#ifndef RENS_FIX_ONLY_LOWER
749        if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
750            value = floor(value + 0.5);
751            if (value == lower || value == upper)
752                numberAtBound++;
753            newSolver->setColLower(iColumn, value);
754            newSolver->setColUpper(iColumn, value);
755            numberFixed++;
756        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
757            numberTightened++;
758            newSolver->setColLower(iColumn, floor(value));
759            newSolver->setColUpper(iColumn, ceil(value));
760        }
761#else
762        if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
763                floor(value + 0.5) == lower &&
764            djValue > djTolerance ) {
765          value = floor(value + 0.5);
766          numberAtBound++;
767          newSolver->setColLower(iColumn, value);
768          newSolver->setColUpper(iColumn, value);
769          numberFixed++;
770        } else if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
771                floor(value + 0.5) == upper &&
772                   -djValue > djTolerance && (djTolerance > 0.0||type==2)) {
773          value = floor(value + 0.5);
774          numberAtBound++;
775          newSolver->setColLower(iColumn, value);
776          newSolver->setColUpper(iColumn, value);
777          numberFixed++;
778        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0 &&
779                   djTolerance <0.0) {
780            numberTightened++;
781            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
782                value = floor(value + 0.5);
783                if (value < upper) {
784                    newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
785                    newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
786                } else {
787                    newSolver->setColLower(iColumn, upper - 1.0);
788                }
789            } else {
790                newSolver->setColLower(iColumn, floor(value));
791                newSolver->setColUpper(iColumn, ceil(value));
792            }
793        }
794#endif
795    }
796    delete [] dj;
797    if (numberFixed > numberIntegers / 5) {
798        if ( numberFixed < numberColumns / 5) {
799#define RENS_FIX_CONTINUOUS
800#ifdef RENS_FIX_CONTINUOUS
801            const double * colLower = newSolver->getColLower();
802            //const double * colUpper = newSolver->getColUpper();
803            int nAtLb = 0;
804            double sumDj = 0.0;
805            const double * dj = newSolver->getReducedCost();
806            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
807                if (!newSolver->isInteger(iColumn)) {
808                    double value = currentSolution[iColumn];
809                    if (value < colLower[iColumn] + 1.0e-8) {
810                        double djValue = dj[iColumn] * direction;
811                        nAtLb++;
812                        sumDj += djValue;
813                    }
814                }
815            }
816            if (nAtLb) {
817                // fix some continuous
818                double * sort = new double[nAtLb];
819                int * which = new int [nAtLb];
820                double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
821                int nFix2 = 0;
822                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
823                    if (!newSolver->isInteger(iColumn)) {
824                        double value = currentSolution[iColumn];
825                        if (value < colLower[iColumn] + 1.0e-8) {
826                            double djValue = dj[iColumn] * direction;
827                            if (djValue > threshold) {
828                                sort[nFix2] = -djValue;
829                                which[nFix2++] = iColumn;
830                            }
831                        }
832                    }
833                }
834                CoinSort_2(sort, sort + nFix2, which);
835                nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
836                for (int i = 0; i < nFix2; i++) {
837                    int iColumn = which[i];
838                    newSolver->setColUpper(iColumn, colLower[iColumn]);
839                }
840                delete [] sort;
841                delete [] which;
842#ifdef CLP_INVESTIGATE2
843                printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
844                       numberFixed, numberTightened, numberAtBound, nFix2);
845#endif
846            }
847#endif
848        }
849#ifdef COIN_DEVELOP
850        printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
851#endif
852        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
853                                         model_->getCutoff(), "CbcHeuristicRENS");
854        if (returnCode < 0) {
855            returnCode = 0; // returned on size
856#ifdef RENS_FIX_CONTINUOUS
857            if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
858                const double * colLower = newSolver->getColLower();
859                //const double * colUpper = newSolver->getColUpper();
860                int nAtLb = 0;
861                double sumDj = 0.0;
862                const double * dj = newSolver->getReducedCost();
863                double direction = newSolver->getObjSense();
864                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
865                    if (!newSolver->isInteger(iColumn)) {
866                        double value = currentSolution[iColumn];
867                        if (value < colLower[iColumn] + 1.0e-8) {
868                            double djValue = dj[iColumn] * direction;
869                            nAtLb++;
870                            sumDj += djValue;
871                        }
872                    }
873                }
874                if (nAtLb) {
875                    // fix some continuous
876                    double * sort = new double[nAtLb];
877                    int * which = new int [nAtLb];
878                    double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
879                    int nFix2 = 0;
880                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
881                        if (!newSolver->isInteger(iColumn)) {
882                            double value = currentSolution[iColumn];
883                            if (value < colLower[iColumn] + 1.0e-8) {
884                                double djValue = dj[iColumn] * direction;
885                                if (djValue > threshold) {
886                                    sort[nFix2] = -djValue;
887                                    which[nFix2++] = iColumn;
888                                }
889                            }
890                        }
891                    }
892                    CoinSort_2(sort, sort + nFix2, which);
893                    nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
894                    for (int i = 0; i < nFix2; i++) {
895                        int iColumn = which[i];
896                        newSolver->setColUpper(iColumn, colLower[iColumn]);
897                    }
898                    delete [] sort;
899                    delete [] which;
900#ifdef CLP_INVESTIGATE2
901                    printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
902                           numberFixed, numberTightened, numberAtBound, nFix2);
903#endif
904                }
905                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
906                                                 model_->getCutoff(), "CbcHeuristicRENS");
907#endif
908            }
909        }
910        //printf("return code %d",returnCode);
911        if ((returnCode&2) != 0) {
912            // could add cut
913            returnCode &= ~2;
914#ifdef COIN_DEVELOP
915            if (!numberTightened && numberFixed == numberAtBound)
916                printf("could add cut with %d elements\n", numberFixed);
917#endif
918        } else {
919            //printf("\n");
920        }
921    }
922    //delete [] whichRow;
923    //delete [] contribution;
924    delete newSolver;
925    fractionSmall_ = saveFractionSmall;
926    return returnCode;
927}
928// update model
929void CbcHeuristicRENS::setModel(CbcModel * model)
930{
931    model_ = model;
932}
933
Note: See TracBrowser for help on using the repository browser.