source: stable/2.8/Cbc/src/CbcHeuristicRENS.cpp @ 1902

Last change on this file since 1902 was 1902, checked in by stefan, 6 years ago

sync with trunk rev 1901

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.8 KB
Line 
1// $Id: CbcHeuristicRENS.cpp 1902 2013-04-10 16:58:16Z stefan $
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            COIN_DETAIL_PRINT(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 || returnCode == 0) {
855#ifdef RENS_FIX_CONTINUOUS
856            if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
857                const double * colLower = newSolver->getColLower();
858                //const double * colUpper = newSolver->getColUpper();
859                int nAtLb = 0;
860                double sumDj = 0.0;
861                const double * dj = newSolver->getReducedCost();
862                double direction = newSolver->getObjSense();
863                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
864                    if (!newSolver->isInteger(iColumn)) {
865                        double value = currentSolution[iColumn];
866                        if (value < colLower[iColumn] + 1.0e-8) {
867                            double djValue = dj[iColumn] * direction;
868                            nAtLb++;
869                            sumDj += djValue;
870                        }
871                    }
872                }
873                if (nAtLb) {
874                    // fix some continuous
875                    double * sort = new double[nAtLb];
876                    int * which = new int [nAtLb];
877                    double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
878                    int nFix2 = 0;
879                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
880                        if (!newSolver->isInteger(iColumn)) {
881                            double value = currentSolution[iColumn];
882                            if (value < colLower[iColumn] + 1.0e-8) {
883                                double djValue = dj[iColumn] * direction;
884                                if (djValue > threshold) {
885                                    sort[nFix2] = -djValue;
886                                    which[nFix2++] = iColumn;
887                                }
888                            }
889                        }
890                    }
891                    CoinSort_2(sort, sort + nFix2, which);
892                    nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
893                    for (int i = 0; i < nFix2; i++) {
894                        int iColumn = which[i];
895                        newSolver->setColUpper(iColumn, colLower[iColumn]);
896                    }
897                    delete [] sort;
898                    delete [] which;
899#ifdef CLP_INVESTIGATE2
900                    printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
901                           numberFixed, numberTightened, numberAtBound, nFix2);
902#endif
903                }
904                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
905                                                 model_->getCutoff(), "CbcHeuristicRENS");
906            }
907#endif
908            if (returnCode < 0 || returnCode == 0) {
909                // Do passes fixing up those >0.9 and
910                // down those < 0.05
911#define RENS_PASS 3
912              //#define KEEP_GOING
913#ifdef KEEP_GOING
914                double * saveLower = CoinCopyOfArray(colLower,numberColumns);
915                double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
916                bool badPass=false;
917                int nSolved=0;
918#endif
919                for (int iPass=0;iPass<RENS_PASS;iPass++) {
920                  int nFixed=0;
921                  int nFixedAlready=0;
922                  int nFixedContinuous=0;
923                  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
924                    if (colUpper[iColumn]>colLower[iColumn]) {
925                      if (newSolver->isInteger(iColumn)) {
926                        double value = currentSolution[iColumn];
927                        double fixTo = floor(value+0.1);
928                        if (fixTo>value || value-fixTo < 0.05) {
929                          // above 0.9 or below 0.05
930                          nFixed++;
931                          newSolver->setColLower(iColumn, fixTo);
932                          newSolver->setColUpper(iColumn, fixTo);
933                        }
934                      }
935                    } else if (newSolver->isInteger(iColumn)) {
936                      nFixedAlready++;
937                    } else {
938                      nFixedContinuous++;
939                    }
940                  }
941#ifdef CLP_INVESTIGATE2
942                  printf("%d more integers fixed (total %d) plus %d continuous\n",
943                         nFixed,nFixed+nFixedAlready,nFixedContinuous);
944#endif
945#ifdef KEEP_GOING
946                  if (nFixed) {
947                    newSolver->resolve();
948                    if (!newSolver->isProvenOptimal()) {
949                      badPass=true;
950                      break;
951                    } else {
952                      nSolved++;
953                      memcpy(saveLower,colLower,numberColumns*sizeof(double));
954                      memcpy(saveUpper,colUpper,numberColumns*sizeof(double));
955                    }
956                  } else {
957                    break;
958                  }
959#else
960                  if (nFixed) {
961                    newSolver->resolve();
962                    if (!newSolver->isProvenOptimal()) {
963                      returnCode=0;
964                      break;
965                    }
966                    returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
967                                                     model_->getCutoff(), "CbcHeuristicRENS");
968                  } else {
969                    returnCode=0;
970                  }
971                  if (returnCode>=0)
972                    break;
973                }
974                if (returnCode < 0) 
975                returnCode = 0; // returned on size
976#endif
977            }
978#ifdef KEEP_GOING
979            if (badPass) {
980              newSolver->setColLower(saveLower);
981              newSolver->setColUpper(saveUpper);
982              newSolver->resolve();
983            }
984            delete [] saveLower;
985            delete [] saveUpper;
986            if (nSolved)
987              returnCode = 
988                smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
989                                    model_->getCutoff(), "CbcHeuristicRENS");
990            else
991              returnCode=0;
992        }
993#endif
994        }
995        //printf("return code %d",returnCode);
996        if ((returnCode&2) != 0) {
997            // could add cut
998            returnCode &= ~2;
999#ifdef COIN_DEVELOP
1000            if (!numberTightened && numberFixed == numberAtBound)
1001                printf("could add cut with %d elements\n", numberFixed);
1002#endif
1003        } else {
1004            //printf("\n");
1005        }
1006    }
1007    //delete [] whichRow;
1008    //delete [] contribution;
1009    delete newSolver;
1010    fractionSmall_ = saveFractionSmall;
1011    return returnCode;
1012}
1013// update model
1014void CbcHeuristicRENS::setModel(CbcModel * model)
1015{
1016    model_ = model;
1017}
1018
Note: See TracBrowser for help on using the repository browser.