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

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

add some more heuristics

File size: 28.7 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 = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
337            char * mark = new char [numberColumns];
338            char * nonzero = new char [numberColumns];
339            double factor=CoinMax(1.000001,fractionSmall_);
340            fractionSmall_ = 0.5;
341            // loosen up
342            for (int i=0;i<numberRows;i++) {
343              if (bestDj[i]>=1.0e30) {
344                newSolver->setRowUpper(i,factor*saveUpper[i]);
345              }
346            }
347            newSolver->resolve();
348            const double * solution = newSolver->getColSolution();
349            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
350              mark[iColumn]=0;
351              nonzero[iColumn]=0;
352              if (colUpper[iColumn]>colLower[iColumn]&&
353                  solution[iColumn]>0.9999)
354                mark[iColumn]=1;
355              else if (solution[iColumn]>0.00001)
356                nonzero[iColumn]=1;
357            }
358            // slightly small
359            for (int i=0;i<numberRows;i++) {
360              if (bestDj[i]>=1.0e30) {
361                newSolver->setRowUpper(i,saveUpper[i]*0.9999);
362              }
363            }
364            newSolver->resolve();
365            int nCheck=2;
366            if (newSolver->isProvenOptimal()) {
367              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
368                if (colUpper[iColumn]>colLower[iColumn]&&
369                    solution[iColumn]>0.9999)
370                  mark[iColumn]++;
371                else if (solution[iColumn]>0.00001)
372                  nonzero[iColumn]=1;
373              }
374            } else {
375              nCheck=1;
376            }
377            // correct values
378            for (int i=0;i<numberRows;i++) {
379              if (bestDj[i]>=1.0e30) {
380                newSolver->setRowUpper(i,saveUpper[i]);
381              }
382            }
383            newSolver->resolve();
384            int nFixed=0;
385            int nFixedToZero=0;
386            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
387              if (colUpper[iColumn]>colLower[iColumn]) {
388                if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) {
389                  newSolver->setColLower(iColumn,1.0);
390                  nFixed++;
391                } else if (!mark[iColumn]&&!nonzero[iColumn]&&
392                           columnLength[iColumn]>1&&solution[iColumn]<0.00001) {
393                  newSolver->setColUpper(iColumn,0.0);
394                  nFixedToZero++;
395                }
396              }
397            }
398            char line[100];
399            sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)",
400                    heuristicName(),
401                    nFixed,nFixedToZero);
402            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
403              << line
404              << CoinMessageEol;
405            delete [] mark;
406            delete []nonzero;
407            delete [] saveUpper;
408            numberFixed=numberColumns;
409            djTolerance = 1.0e30;
410          } else if (type==11) {
411            double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
412            char * mark = new char [numberColumns];
413            char * nonzero = new char [numberColumns];
414            // save basis and solution
415            CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(newSolver->getWarmStart()) ;
416            assert(basis != NULL);
417            double * saveSolution = 
418              CoinCopyOfArray(newSolver->getColSolution(), 
419                              numberColumns);
420            double factors[] = {1.1,1.05,1.01,0.98};
421            int nPass = (sizeof(factors)/sizeof(double))-1;
422            double factor=factors[0];
423            double proportion = fractionSmall_;
424            fractionSmall_ = 0.5;
425            // loosen up
426            for (int i=0;i<numberRows;i++) {
427              if (bestDj[i]>=1.0e30) {
428                newSolver->setRowUpper(i,factor*saveUpper[i]);
429              }
430            }
431            bool takeHint;
432            OsiHintStrength strength;
433            newSolver->getHintParam(OsiDoDualInResolve, takeHint, strength);
434            newSolver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
435            newSolver->resolve();
436            newSolver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
437            const double * solution = newSolver->getColSolution();
438            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
439              mark[iColumn]=0;
440              nonzero[iColumn]=0;
441              if (colUpper[iColumn]>colLower[iColumn]&&
442                  solution[iColumn]>0.9999)
443                mark[iColumn]=1;
444              else if (solution[iColumn]>0.00001)
445                nonzero[iColumn]=1;
446            }
447            int nCheck=2;
448            for (int iPass=0;iPass<nPass;iPass++) {
449              // smaller
450              factor = factors[iPass+1];
451              for (int i=0;i<numberRows;i++) {
452                if (bestDj[i]>=1.0e30) {
453                  newSolver->setRowUpper(i,saveUpper[i]*factor);
454                }
455              }
456              newSolver->resolve();
457              if (newSolver->isProvenOptimal()) {
458                nCheck++;
459                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
460                  if (colUpper[iColumn]>colLower[iColumn]&&
461                      solution[iColumn]>0.9999)
462                    mark[iColumn]++;
463                  else if (solution[iColumn]>0.00001)
464                    nonzero[iColumn]++;
465                }
466              }
467            }
468            // correct values
469            for (int i=0;i<numberRows;i++) {
470              if (bestDj[i]>=1.0e30) {
471                newSolver->setRowUpper(i,saveUpper[i]);
472              }
473            }
474            newSolver->setColSolution(saveSolution);
475            delete [] saveSolution;
476            newSolver->setWarmStart(basis);
477            delete basis ;
478            newSolver->setHintParam(OsiDoDualInResolve, takeHint, strength);
479            newSolver->resolve();
480            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
481              if (colUpper[iColumn]>colLower[iColumn]&&
482                  solution[iColumn]>0.9999)
483                mark[iColumn]++;
484              else if (solution[iColumn]>0.00001)
485                nonzero[iColumn]++;
486            }
487            int nFixed=0;
488            int numberSetsToFix = static_cast<int>(nSOS*(1.0-proportion));
489            int * mixed = new int[numberRows];
490            memset(mixed,0,numberRows*sizeof(int));
491            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
492              if (colUpper[iColumn]>colLower[iColumn]) {
493                int iSOS=-1;
494                for (CoinBigIndex j = columnStart[iColumn];
495                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
496                  int iRow = row[j];
497                  if (bestDj[iRow]<1.0e25) {
498                    iSOS=iRow;
499                    break;
500                  }
501                }
502                if (iSOS>=0) {
503                  int numberTimesAtOne = mark[iColumn];
504                  int numberTimesNonZero = nonzero[iColumn]+
505                    numberTimesAtOne;
506                  if (numberTimesAtOne<nCheck&&
507                      numberTimesNonZero) {
508                    mixed[iSOS]+=
509                      CoinMin(numberTimesNonZero,
510                              nCheck-numberTimesNonZero);
511                  } 
512                }
513              }
514            }
515            int nFix=0;
516            for (int i=0;i<numberRows;i++) {
517              if (bestDj[i]<1.0e25) {
518                sort[nFix] = -bestDj[i]+1.0e8*mixed[i];
519                mixed[nFix++]=i;
520              }
521            }
522            CoinSort_2(sort,sort+nFix,mixed);
523            nFix = CoinMin(nFix,numberSetsToFix);
524            memset(sort,0,sizeof(double)*numberRows);
525            for (int i=0;i<nFix;i++)
526              sort[mixed[i]]=1.0;
527            delete [] mixed;
528            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
529              if (colUpper[iColumn]>colLower[iColumn]) {
530                if (solution[iColumn]>0.9999) {
531                  int iSOS=-1;
532                  for (CoinBigIndex j = columnStart[iColumn];
533                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
534                    int iRow = row[j];
535                    if (bestDj[iRow]<1.0e25) {
536                      iSOS=iRow;
537                      break;
538                    }
539                  }
540                  if (iSOS>=0&&sort[iSOS]) {
541                    newSolver->setColLower(iColumn,1.0);
542                    nFixed++;
543                  }
544                }
545              }
546            }
547            char line[100];
548            sprintf(line,"Heuristic %s fixed %d to one",
549                    heuristicName(),
550                    nFixed);
551            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
552              << line
553              << CoinMessageEol;
554            delete [] mark;
555            delete [] nonzero;
556            delete [] saveUpper;
557            numberFixed=numberColumns;
558            djTolerance = 1.0e30;
559          }
560        }
561        delete basis;
562        delete [] sort;
563        delete [] bestDj;
564        if (10*nSOS<=8*numberRows) {
565          // give up
566          delete [] contribution;
567          delete newSolver;
568          return 0;
569        }
570      }
571    }
572    // Do dj to get right number
573    if (type==4||type==6||(type>7&&type<10)) {
574      double * sort = new double [numberColumns];
575      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
576        sort[iColumn]=1.0e30;
577        if (colUpper[iColumn]>colLower[iColumn]) {
578          sort[iColumn] = fabs(dj[iColumn]);
579        }
580      }
581      std::sort(sort,sort+numberColumns);
582      int last = static_cast<int>(numberColumns*fractionSmall_);
583      djTolerance = CoinMax(sort[last],1.0e-5);
584      delete [] sort;
585    } else if (type==12) {
586      // Do layered in a different way
587      int numberRows = solver->getNumRows();
588      // Column copy
589      const CoinPackedMatrix * matrix = newSolver->getMatrixByCol();
590      const double * element = matrix->getElements();
591      const int * row = matrix->getIndices();
592      const CoinBigIndex * columnStart = matrix->getVectorStarts();
593      const int * columnLength = matrix->getVectorLengths();
594      int * whichRow = new int[numberRows];
595      int * whichSet = new int [numberColumns];
596      int nSOS=0;
597      for (int i=0;i<numberRows;i++) {
598        whichRow[i]=0;
599        if (!contribution[i])
600          nSOS++;
601      }
602      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
603        whichSet[iColumn]=-2;
604        if (colUpper[iColumn]>colLower[iColumn]) {
605          CoinBigIndex j;
606          double sum=0.0;
607          int iSOS=-1;
608          int n=0;
609          for (j = columnStart[iColumn];
610               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
611            int iRow = row[j];
612            if (contribution[iRow]>=0.0) {
613              iSOS=iRow;
614              n++;
615            } else {
616              sum += fabs(element[j]);
617            }
618          }
619          if (n>1)
620            printf("Too many SOS entries (%d) for column %d\n",
621                   n,iColumn);
622          if (sum) {
623            assert (iSOS>=0);
624            contribution[iSOS] += sum;
625            whichRow[iSOS]++;
626            whichSet[iColumn]=iSOS;
627          } else {
628            whichSet[iColumn]=iSOS+numberRows;
629          }
630        }
631      }
632      int * chunk = new int [numberRows];
633      for (int i=0;i<numberRows;i++) {
634        chunk[i]=-1;
635        if (whichRow[i]) {
636          contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]);
637        } else {
638          contribution[i] = COIN_DBL_MAX;
639        }
640        whichRow[i]=i;
641      }
642      newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
643      double * saveLower = CoinCopyOfArray(colLower,numberColumns);
644      double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
645      CoinSort_2(contribution,contribution+numberRows,whichRow);
646      // Set do nothing solution
647      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
648        if(whichSet[iColumn]>=numberRows)
649          newSolver->setColLower(iColumn,1.0);
650      }
651      newSolver->resolve();
652      int nChunk = (nSOS+9)/10;
653      int nPass=0;
654      int inChunk=0;
655      for (int i=0;i<nSOS;i++) {
656        chunk[whichRow[i]]=nPass;
657        inChunk++;
658        if (inChunk==nChunk) {
659          inChunk=0;
660          // last two together
661          if (i+nChunk<nSOS)
662            nPass++;
663        }
664      }
665      // adjust
666      nPass++;
667      for (int iPass=0;iPass<nPass;iPass++) {
668        // fix last chunk and unfix this chunk
669        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
670          int iSOS = whichSet[iColumn];
671          if (iSOS>=0) {
672            if (iSOS>=numberRows)
673              iSOS-=numberRows;
674            if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) {
675              newSolver->setColLower(iColumn,1.0);
676            } else if (chunk[iSOS]==iPass) {
677              newSolver->setColLower(iColumn,saveLower[iColumn]);
678              newSolver->setColUpper(iColumn,saveUpper[iColumn]);
679            }
680          }
681        }
682        // solve
683        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
684                                         model_->getCutoff(), "CbcHeuristicRENS");
685        if (returnCode < 0) {
686            returnCode = 0; // returned on size
687            break;
688        } else if ((returnCode&1)==0) {
689          // no good
690          break;
691        }
692      }
693      if ((returnCode&2) != 0) {
694        // could add cut
695        returnCode &= ~2;
696      }
697      delete [] chunk;
698      delete [] saveLower;
699      delete [] saveUpper;
700      delete [] whichRow;
701      delete [] whichSet;
702      delete [] contribution;
703      delete newSolver;
704      return returnCode;
705    }
706   
707    double primalTolerance;
708    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
709
710    int i;
711    int numberTightened = 0;
712    int numberAtBound = 0;
713    int numberContinuous = numberColumns - numberIntegers;
714
715    for (i = 0; i < numberIntegers; i++) {
716        int iColumn = integerVariable[i];
717        double value = currentSolution[iColumn];
718        double lower = colLower[iColumn];
719        double upper = colUpper[iColumn];
720        value = CoinMax(value, lower);
721        value = CoinMin(value, upper);
722        double djValue=dj[iColumn]*direction;
723#define RENS_FIX_ONLY_LOWER
724#ifndef RENS_FIX_ONLY_LOWER
725        if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
726            value = floor(value + 0.5);
727            if (value == lower || value == upper)
728                numberAtBound++;
729            newSolver->setColLower(iColumn, value);
730            newSolver->setColUpper(iColumn, value);
731            numberFixed++;
732        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
733            numberTightened++;
734            newSolver->setColLower(iColumn, floor(value));
735            newSolver->setColUpper(iColumn, ceil(value));
736        }
737#else
738        if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
739                floor(value + 0.5) == lower &&
740            djValue > djTolerance ) {
741          value = floor(value + 0.5);
742          numberAtBound++;
743          newSolver->setColLower(iColumn, value);
744          newSolver->setColUpper(iColumn, value);
745          numberFixed++;
746        } else if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
747                floor(value + 0.5) == upper &&
748                   -djValue > djTolerance && (djTolerance > 0.0||type==2)) {
749          value = floor(value + 0.5);
750          numberAtBound++;
751          newSolver->setColLower(iColumn, value);
752          newSolver->setColUpper(iColumn, value);
753          numberFixed++;
754        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0 &&
755                   djTolerance <0.0) {
756            numberTightened++;
757            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
758                value = floor(value + 0.5);
759                if (value < upper) {
760                    newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
761                    newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
762                } else {
763                    newSolver->setColLower(iColumn, upper - 1.0);
764                }
765            } else {
766                newSolver->setColLower(iColumn, floor(value));
767                newSolver->setColUpper(iColumn, ceil(value));
768            }
769        }
770#endif
771    }
772    delete [] dj;
773    if (numberFixed > numberIntegers / 5) {
774        if ( numberFixed < numberColumns / 5) {
775#define RENS_FIX_CONTINUOUS
776#ifdef RENS_FIX_CONTINUOUS
777            const double * colLower = newSolver->getColLower();
778            //const double * colUpper = newSolver->getColUpper();
779            int nAtLb = 0;
780            double sumDj = 0.0;
781            const double * dj = newSolver->getReducedCost();
782            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
783                if (!newSolver->isInteger(iColumn)) {
784                    double value = currentSolution[iColumn];
785                    if (value < colLower[iColumn] + 1.0e-8) {
786                        double djValue = dj[iColumn] * direction;
787                        nAtLb++;
788                        sumDj += djValue;
789                    }
790                }
791            }
792            if (nAtLb) {
793                // fix some continuous
794                double * sort = new double[nAtLb];
795                int * which = new int [nAtLb];
796                double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
797                int nFix2 = 0;
798                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
799                    if (!newSolver->isInteger(iColumn)) {
800                        double value = currentSolution[iColumn];
801                        if (value < colLower[iColumn] + 1.0e-8) {
802                            double djValue = dj[iColumn] * direction;
803                            if (djValue > threshold) {
804                                sort[nFix2] = -djValue;
805                                which[nFix2++] = iColumn;
806                            }
807                        }
808                    }
809                }
810                CoinSort_2(sort, sort + nFix2, which);
811                nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
812                for (int i = 0; i < nFix2; i++) {
813                    int iColumn = which[i];
814                    newSolver->setColUpper(iColumn, colLower[iColumn]);
815                }
816                delete [] sort;
817                delete [] which;
818#ifdef CLP_INVESTIGATE2
819                printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
820                       numberFixed, numberTightened, numberAtBound, nFix2);
821#endif
822            }
823#endif
824        }
825#ifdef COIN_DEVELOP
826        printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
827#endif
828        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
829                                         model_->getCutoff(), "CbcHeuristicRENS");
830        if (returnCode < 0) {
831            returnCode = 0; // returned on size
832#ifdef RENS_FIX_CONTINUOUS
833            if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
834                const double * colLower = newSolver->getColLower();
835                //const double * colUpper = newSolver->getColUpper();
836                int nAtLb = 0;
837                double sumDj = 0.0;
838                const double * dj = newSolver->getReducedCost();
839                double direction = newSolver->getObjSense();
840                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
841                    if (!newSolver->isInteger(iColumn)) {
842                        double value = currentSolution[iColumn];
843                        if (value < colLower[iColumn] + 1.0e-8) {
844                            double djValue = dj[iColumn] * direction;
845                            nAtLb++;
846                            sumDj += djValue;
847                        }
848                    }
849                }
850                if (nAtLb) {
851                    // fix some continuous
852                    double * sort = new double[nAtLb];
853                    int * which = new int [nAtLb];
854                    double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
855                    int nFix2 = 0;
856                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
857                        if (!newSolver->isInteger(iColumn)) {
858                            double value = currentSolution[iColumn];
859                            if (value < colLower[iColumn] + 1.0e-8) {
860                                double djValue = dj[iColumn] * direction;
861                                if (djValue > threshold) {
862                                    sort[nFix2] = -djValue;
863                                    which[nFix2++] = iColumn;
864                                }
865                            }
866                        }
867                    }
868                    CoinSort_2(sort, sort + nFix2, which);
869                    nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
870                    for (int i = 0; i < nFix2; i++) {
871                        int iColumn = which[i];
872                        newSolver->setColUpper(iColumn, colLower[iColumn]);
873                    }
874                    delete [] sort;
875                    delete [] which;
876#ifdef CLP_INVESTIGATE2
877                    printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
878                           numberFixed, numberTightened, numberAtBound, nFix2);
879#endif
880                }
881                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
882                                                 model_->getCutoff(), "CbcHeuristicRENS");
883#endif
884            }
885        }
886        //printf("return code %d",returnCode);
887        if ((returnCode&2) != 0) {
888            // could add cut
889            returnCode &= ~2;
890#ifdef COIN_DEVELOP
891            if (!numberTightened && numberFixed == numberAtBound)
892                printf("could add cut with %d elements\n", numberFixed);
893#endif
894        } else {
895            //printf("\n");
896        }
897    }
898    //delete [] whichRow;
899    //delete [] contribution;
900    delete newSolver;
901    fractionSmall_ = saveFractionSmall;
902    return returnCode;
903}
904// update model
905void CbcHeuristicRENS::setModel(CbcModel * model)
906{
907    model_ = model;
908}
909
Note: See TracBrowser for help on using the repository browser.