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

Last change on this file since 2105 was 2105, checked in by forrest, 3 years ago

mostly reporting options plus a few tweaks

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.6 KB
Line 
1// $Id: CbcHeuristicRENS.cpp 2105 2015-01-05 13:11:11Z forrest $
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// Constructor with model - assumed before cuts
36
37CbcHeuristicRENS::CbcHeuristicRENS(CbcModel & model)
38        : CbcHeuristic(model)
39{
40    numberTries_ = 0;
41    rensType_ = 0;
42    whereFrom_ = 256 + 1;
43}
44
45// Destructor
46CbcHeuristicRENS::~CbcHeuristicRENS ()
47{
48}
49
50// Clone
51CbcHeuristic *
52CbcHeuristicRENS::clone() const
53{
54    return new CbcHeuristicRENS(*this);
55}
56
57// Assignment operator
58CbcHeuristicRENS &
59CbcHeuristicRENS::operator=( const CbcHeuristicRENS & rhs)
60{
61    if (this != &rhs) {
62        CbcHeuristic::operator=(rhs);
63        numberTries_ = rhs.numberTries_;
64        rensType_ = rhs.rensType_;
65    }
66    return *this;
67}
68
69// Copy constructor
70CbcHeuristicRENS::CbcHeuristicRENS(const CbcHeuristicRENS & rhs)
71        :
72        CbcHeuristic(rhs),
73        numberTries_(rhs.numberTries_),
74        rensType_(rhs.rensType_)
75{
76}
77// Resets stuff if model changes
78void
79CbcHeuristicRENS::resetModel(CbcModel * )
80{
81}
82int
83CbcHeuristicRENS::solution(double & solutionValue,
84                           double * betterSolution)
85{
86    int returnCode = 0;
87    const double * bestSolution = model_->bestSolution();
88    bool returnNow=false;
89    if ((numberTries_&&(rensType_&16)==0) || numberTries_>1 || (when() < 2 && bestSolution))
90      returnNow=true;;
91    // If 32 bit set then do once with bestSolution
92    if ((rensType_&32)!=0&&bestSolution)
93      returnNow=false;
94    if (returnNow)
95        return 0;
96    // switch off next time if bestSolution
97    if (bestSolution) {
98      if ((rensType_&32)!=0)
99        rensType_ &= ~32; // switch off but leave bestSolution
100      else
101        bestSolution=NULL; // null bestSolution so won't use
102    }
103    numberTries_++;
104#ifdef HEURISTIC_INFORM
105    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
106           heuristicName(),numRuns_,numCouldRun_,when_);
107#endif
108    double saveFractionSmall=fractionSmall_;
109    OsiSolverInterface * solver = model_->solver();
110
111    int numberIntegers = model_->numberIntegers();
112    const int * integerVariable = model_->integerVariable();
113
114    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
115    const double * currentSolution = newSolver->getColSolution();
116    int type = rensType_&15;
117    if (type<12)
118      newSolver->resolve();
119    double direction = newSolver->getObjSense();
120    double cutoff=model_->getCutoff();
121    newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
122    //cutoff *= direction;
123    double gap = cutoff - newSolver->getObjValue() * direction ;
124    double tolerance;
125    newSolver->getDblParam(OsiDualTolerance, tolerance) ;
126    if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<12) {
127      gap += 100.0 * tolerance;
128      int nFix = newSolver->reducedCostFix(gap);
129      if (nFix) {
130        char line [200];
131        sprintf(line, "Reduced cost fixing fixed %d variables", nFix);
132        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
133          << line
134          << CoinMessageEol;
135      }
136    } else if (type<12) {
137      return 0; // finished?
138    }
139    int numberColumns = solver->getNumCols();
140    double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
141    double djTolerance = (type!=1) ? -1.0e30 : 1.0e-4;
142    const double * colLower = newSolver->getColLower();
143    const double * colUpper = newSolver->getColUpper();
144    double * contribution = NULL;
145    int numberFixed = 0;
146    if (type==3) {
147      double total=0.0;
148      int n=0;
149      CoinWarmStartBasis * basis =
150        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
151      if (basis&&basis->getNumArtificial()) {
152        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
153          if (colUpper[iColumn]>colLower[iColumn]&&
154              basis->getStructStatus(iColumn) !=
155              CoinWarmStartBasis::basic) {
156            n++;
157            total += fabs(dj[iColumn]);
158          }
159        }
160        if (n)
161          djTolerance = (0.01*total)/static_cast<double>(n);
162        delete basis;
163      }
164    } else if (type>=5&&type<=12) {
165      /* 5 fix sets at one
166         6 fix on dj but leave unfixed SOS slacks
167         7 fix sets at one but use pi
168         8 fix all at zero but leave unfixed SOS slacks
169         9 as 8 but only fix all at zero if just one in set nonzero
170         10 fix all "stable" ones
171         11 fix all "stable" ones - approach 2
172         12 layered approach
173      */
174      // SOS type fixing
175      bool fixSets = (type==5)||(type==7)||(type==10)||(type==11);
176      CoinWarmStartBasis * basis =
177        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
178      if (basis&&basis->getNumArtificial()) {
179        //const double * rowLower = solver->getRowLower();
180        const double * rowUpper = solver->getRowUpper();
181       
182        int numberRows = solver->getNumRows();
183        // Column copy
184        const CoinPackedMatrix * matrix = solver->getMatrixByCol();
185        const double * element = matrix->getElements();
186        const int * row = matrix->getIndices();
187        const CoinBigIndex * columnStart = matrix->getVectorStarts();
188        const int * columnLength = matrix->getVectorLengths();
189        double * bestDj = new double [numberRows];
190        for (int i=0;i<numberRows;i++) {
191          if (rowUpper[i]==1.0)
192            bestDj[i]=1.0e20;
193          else
194            bestDj[i]=1.0e30;
195        }
196        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
197          if (colUpper[iColumn]>colLower[iColumn]) {
198            CoinBigIndex j;
199            if (currentSolution[iColumn]>1.0e-6&&
200                currentSolution[iColumn]<0.999999) {
201              for (j = columnStart[iColumn];
202                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
203                int iRow = row[j];
204                if (bestDj[iRow]<1.0e30) {
205                  if (element[j] != 1.0)
206                    bestDj[iRow]=1.0e30;
207                  else
208                    bestDj[iRow]=1.0e25;
209                }
210              }
211            } else if ( basis->getStructStatus(iColumn) !=
212              CoinWarmStartBasis::basic) {
213              for (j = columnStart[iColumn];
214                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
215                int iRow = row[j];
216                if (bestDj[iRow]<1.0e25) {
217                  if (element[j] != 1.0)
218                    bestDj[iRow]=1.0e30;
219                  else
220                    bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]);
221                }
222              }
223            }
224          }
225        }
226        // Just leave one slack in each set
227        {
228          const double * objective = newSolver->getObjCoefficients();
229          int * best = new int [numberRows];
230          double * cheapest = new double[numberRows];
231          for (int i=0;i<numberRows;i++) {
232            best[i]=-1;
233            cheapest[i]=COIN_DBL_MAX;
234          }
235          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
236            if (colUpper[iColumn]>colLower[iColumn]) {
237              if (columnLength[iColumn]==1) {
238                CoinBigIndex j = columnStart[iColumn];
239                int iRow = row[j];
240                if (bestDj[iRow]<1.0e30) {
241                  double obj = direction*objective[iColumn];
242                  if (obj<cheapest[iRow]) {
243                    cheapest[iRow]=obj;
244                    best[iRow]=iColumn;
245                  }
246                }
247              }
248            }
249          }
250          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
251            if (colUpper[iColumn]>colLower[iColumn]) {
252              if (columnLength[iColumn]==1) {
253                CoinBigIndex j = columnStart[iColumn];
254                int iRow = row[j];
255                if (bestDj[iRow]<1.0e30) {
256                  if (best[iRow]!=-1&&iColumn!=best[iRow]) {
257                    newSolver->setColUpper(iColumn,0.0);
258                  }
259                }
260              }
261            }
262          }
263          delete [] best;
264          delete [] cheapest;
265        }
266        int nSOS=0;
267        double * sort = new double [numberRows];
268        const double * pi = newSolver->getRowPrice();
269        if (type==12) {
270          contribution = new double [numberRows];
271          for (int i=0;i<numberRows;i++) {
272            if (bestDj[i]<1.0e30) 
273              contribution[i]=0.0;
274            else
275              contribution[i]=-1.0;
276          }
277        }
278        for (int i=0;i<numberRows;i++) {
279          if (bestDj[i]<1.0e30) {
280            if (type==5)
281              sort[nSOS++]=bestDj[i];
282            else if (type==7)
283              sort[nSOS++]=-fabs(pi[i]);
284            else
285              sort[nSOS++]=fabs(pi[i]);
286          }
287        }
288        if (10*nSOS>8*numberRows) {
289          if (type<10) {
290            std::sort(sort,sort+nSOS);
291            int last = static_cast<int>(nSOS*0.9*fractionSmall_);
292            double tolerance = sort[last];
293            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
294              if (colUpper[iColumn]>colLower[iColumn]) {
295                CoinBigIndex j;
296                if (currentSolution[iColumn]<=1.0e-6||
297                    currentSolution[iColumn]>=0.999999) {
298                  if (fixSets) {
299                    for (j = columnStart[iColumn];
300                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
301                      int iRow = row[j];
302                      double useDj;
303                      if (type==5) 
304                        useDj = bestDj[iRow];
305                      else if (type==7)
306                        useDj= -fabs(pi[iRow]);
307                      else
308                        useDj= fabs(pi[iRow]);
309                      if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
310                        numberFixed++;
311                        if (currentSolution[iColumn]<=1.0e-6)
312                          newSolver->setColUpper(iColumn,0.0);
313                        else if (currentSolution[iColumn]>=0.999999) 
314                          newSolver->setColLower(iColumn,1.0);
315                      }
316                    }
317                  } else if (columnLength[iColumn]==1) {
318                    // leave more slacks
319                    int iRow = row[columnStart[iColumn]];
320                    if (bestDj[iRow]<1.0e30) {
321                      // fake dj
322                      dj[iColumn] *= 0.000001;
323                    }
324                  } else if (type==8||type==9) {
325                    if (currentSolution[iColumn]<=1.0e-6) {
326                      if (type==8) {
327                        dj[iColumn] *= 1.0e6;
328                      } else {
329                        bool fix=false;
330                        for (j = columnStart[iColumn];
331                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
332                          int iRow = row[j];
333                          if (bestDj[iRow]<1.0e25) {
334                            fix=true;
335                            break;
336                          }
337                        }
338                        if (fix) {
339                          dj[iColumn] *= 1.0e6;
340                        }
341                      }
342                    } else {
343                      dj[iColumn] *= 0.000001;
344                    }
345                  }
346                }
347              }
348            }
349            if (fixSets)
350              djTolerance = 1.0e30;
351          } else if (type==10) {
352            double * saveUpper = new double [numberRows];
353            memset(saveUpper,0,numberRows*sizeof(double));
354            char * mark = new char [numberColumns];
355            char * nonzero = new char [numberColumns];
356            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
357              if (colUpper[iColumn]>colLower[iColumn]) {
358                CoinBigIndex j;
359                for (j = columnStart[iColumn];
360                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
361                  int iRow = row[j];
362                  saveUpper[iRow] += element[j];
363                }
364              }
365            }
366            double sum=0.0;
367            double sumRhs=0.0;
368            const double * rowUpper = newSolver->getRowUpper();
369            for (int i=0;i<numberRows;i++) {
370              if (bestDj[i]>=1.0e30) {
371                sum += saveUpper[i];
372                sumRhs += rowUpper[i];
373              }
374            }
375            double averagePerSet = sum/static_cast<double>(numberRows);
376            // allow this extra
377            double factor = averagePerSet*fractionSmall_*numberRows;
378            factor = 1.0+factor/sumRhs;
379            fractionSmall_ = 0.5;
380            memcpy(saveUpper,rowUpper,numberRows*sizeof(double));
381            // loosen up
382            for (int i=0;i<numberRows;i++) {
383              if (bestDj[i]>=1.0e30) {
384                newSolver->setRowUpper(i,factor*saveUpper[i]);
385              }
386            }
387            newSolver->resolve();
388            const double * solution = newSolver->getColSolution();
389            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
390              mark[iColumn]=0;
391              nonzero[iColumn]=0;
392              if (colUpper[iColumn]>colLower[iColumn]&&
393                  solution[iColumn]>0.9999)
394                mark[iColumn]=1;
395              else if (solution[iColumn]>0.00001)
396                nonzero[iColumn]=1;
397            }
398            // slightly small
399            for (int i=0;i<numberRows;i++) {
400              if (bestDj[i]>=1.0e30) {
401                newSolver->setRowUpper(i,saveUpper[i]*0.9999);
402              }
403            }
404            newSolver->resolve();
405            int nCheck=2;
406            if (newSolver->isProvenOptimal()) {
407              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
408                if (colUpper[iColumn]>colLower[iColumn]&&
409                    solution[iColumn]>0.9999)
410                  mark[iColumn]++;
411                else if (solution[iColumn]>0.00001)
412                  nonzero[iColumn]=1;
413              }
414            } else {
415              nCheck=1;
416            }
417            // correct values
418            for (int i=0;i<numberRows;i++) {
419              if (bestDj[i]>=1.0e30) {
420                newSolver->setRowUpper(i,saveUpper[i]);
421              }
422            }
423            newSolver->resolve();
424            int nFixed=0;
425            int nFixedToZero=0;
426            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
427              if (colUpper[iColumn]>colLower[iColumn]) {
428                if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) {
429                  newSolver->setColLower(iColumn,1.0);
430                  nFixed++;
431                } else if (!mark[iColumn]&&!nonzero[iColumn]&&
432                           columnLength[iColumn]>1&&solution[iColumn]<0.00001) {
433                  newSolver->setColUpper(iColumn,0.0);
434                  nFixedToZero++;
435                }
436              }
437            }
438            char line[100];
439            sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)",
440                    heuristicName(),
441                    nFixed,nFixedToZero);
442            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
443              << line
444              << CoinMessageEol;
445            delete [] mark;
446            delete []nonzero;
447            delete [] saveUpper;
448            numberFixed=numberColumns;
449            djTolerance = 1.0e30;
450          } else if (type==11) {
451            double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
452            char * mark = new char [numberColumns];
453            char * nonzero = new char [numberColumns];
454            // save basis and solution
455            CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(newSolver->getWarmStart()) ;
456            assert(basis != NULL);
457            double * saveSolution = 
458              CoinCopyOfArray(newSolver->getColSolution(), 
459                              numberColumns);
460            double factors[] = {1.1,1.05,1.01,0.98};
461            int nPass = (sizeof(factors)/sizeof(double))-1;
462            double factor=factors[0];
463            double proportion = fractionSmall_;
464            fractionSmall_ = 0.5;
465            // loosen up
466            for (int i=0;i<numberRows;i++) {
467              if (bestDj[i]>=1.0e30) {
468                newSolver->setRowUpper(i,factor*saveUpper[i]);
469              }
470            }
471            bool takeHint;
472            OsiHintStrength strength;
473            newSolver->getHintParam(OsiDoDualInResolve, takeHint, strength);
474            newSolver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
475            newSolver->resolve();
476            newSolver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
477            const double * solution = newSolver->getColSolution();
478            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
479              mark[iColumn]=0;
480              nonzero[iColumn]=0;
481              if (colUpper[iColumn]>colLower[iColumn]&&
482                  solution[iColumn]>0.9999)
483                mark[iColumn]=1;
484              else if (solution[iColumn]>0.00001)
485                nonzero[iColumn]=1;
486            }
487            int nCheck=2;
488            for (int iPass=0;iPass<nPass;iPass++) {
489              // smaller
490              factor = factors[iPass+1];
491              for (int i=0;i<numberRows;i++) {
492                if (bestDj[i]>=1.0e30) {
493                  newSolver->setRowUpper(i,saveUpper[i]*factor);
494                }
495              }
496              newSolver->resolve();
497              if (newSolver->isProvenOptimal()) {
498                nCheck++;
499                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
500                  if (colUpper[iColumn]>colLower[iColumn]&&
501                      solution[iColumn]>0.9999)
502                    mark[iColumn]++;
503                  else if (solution[iColumn]>0.00001)
504                    nonzero[iColumn]++;
505                }
506              }
507            }
508            // correct values
509            for (int i=0;i<numberRows;i++) {
510              if (bestDj[i]>=1.0e30) {
511                newSolver->setRowUpper(i,saveUpper[i]);
512              }
513            }
514            newSolver->setColSolution(saveSolution);
515            delete [] saveSolution;
516            newSolver->setWarmStart(basis);
517            delete basis ;
518            newSolver->setHintParam(OsiDoDualInResolve, takeHint, strength);
519            newSolver->resolve();
520            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
521              if (colUpper[iColumn]>colLower[iColumn]&&
522                  solution[iColumn]>0.9999)
523                mark[iColumn]++;
524              else if (solution[iColumn]>0.00001)
525                nonzero[iColumn]++;
526            }
527            int nFixed=0;
528            int numberSetsToFix = static_cast<int>(nSOS*(1.0-proportion));
529            int * mixed = new int[numberRows];
530            memset(mixed,0,numberRows*sizeof(int));
531            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
532              if (colUpper[iColumn]>colLower[iColumn]) {
533                int iSOS=-1;
534                for (CoinBigIndex j = columnStart[iColumn];
535                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
536                  int iRow = row[j];
537                  if (bestDj[iRow]<1.0e25) {
538                    iSOS=iRow;
539                    break;
540                  }
541                }
542                if (iSOS>=0) {
543                  int numberTimesAtOne = mark[iColumn];
544                  int numberTimesNonZero = nonzero[iColumn]+
545                    numberTimesAtOne;
546                  if (numberTimesAtOne<nCheck&&
547                      numberTimesNonZero) {
548                    mixed[iSOS]+=
549                      CoinMin(numberTimesNonZero,
550                              nCheck-numberTimesNonZero);
551                  } 
552                }
553              }
554            }
555            int nFix=0;
556            for (int i=0;i<numberRows;i++) {
557              if (bestDj[i]<1.0e25) {
558                sort[nFix] = -bestDj[i]+1.0e8*mixed[i];
559                mixed[nFix++]=i;
560              }
561            }
562            CoinSort_2(sort,sort+nFix,mixed);
563            nFix = CoinMin(nFix,numberSetsToFix);
564            memset(sort,0,sizeof(double)*numberRows);
565            for (int i=0;i<nFix;i++)
566              sort[mixed[i]]=1.0;
567            delete [] mixed;
568            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
569              if (colUpper[iColumn]>colLower[iColumn]) {
570                if (solution[iColumn]>0.9999) {
571                  int iSOS=-1;
572                  for (CoinBigIndex j = columnStart[iColumn];
573                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
574                    int iRow = row[j];
575                    if (bestDj[iRow]<1.0e25) {
576                      iSOS=iRow;
577                      break;
578                    }
579                  }
580                  if (iSOS>=0&&sort[iSOS]) {
581                    newSolver->setColLower(iColumn,1.0);
582                    nFixed++;
583                  }
584                }
585              }
586            }
587            char line[100];
588            sprintf(line,"Heuristic %s fixed %d to one (%d sets)",
589                    heuristicName(),
590                    nFixed,nSOS);
591            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
592              << line
593              << CoinMessageEol;
594            delete [] mark;
595            delete [] nonzero;
596            delete [] saveUpper;
597            numberFixed=numberColumns;
598            djTolerance = 1.0e30;
599          }
600        }
601        delete basis;
602        delete [] sort;
603        delete [] bestDj;
604        if (10*nSOS<=8*numberRows) {
605          // give up
606          delete [] contribution;
607          delete newSolver;
608          return 0;
609        }
610      }
611    }
612    // Do dj to get right number
613    if (type==4||type==6||(type>7&&type<10)) {
614      double * sort = new double [numberColumns];
615      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
616        sort[iColumn]=1.0e30;
617        if (colUpper[iColumn]>colLower[iColumn]) {
618          sort[iColumn] = fabs(dj[iColumn]);
619        }
620      }
621      std::sort(sort,sort+numberColumns);
622      int last = static_cast<int>(numberColumns*fractionSmall_);
623      djTolerance = CoinMax(sort[last],1.0e-5);
624      delete [] sort;
625    } else if (type==12) {
626      // Do layered in a different way
627      int numberRows = solver->getNumRows();
628      // Column copy
629      const CoinPackedMatrix * matrix = newSolver->getMatrixByCol();
630      const double * element = matrix->getElements();
631      const int * row = matrix->getIndices();
632      const CoinBigIndex * columnStart = matrix->getVectorStarts();
633      const int * columnLength = matrix->getVectorLengths();
634      int * whichRow = new int[numberRows];
635      int * whichSet = new int [numberColumns];
636      int nSOS=0;
637      for (int i=0;i<numberRows;i++) {
638        whichRow[i]=0;
639        if (!contribution[i])
640          nSOS++;
641      }
642      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
643        whichSet[iColumn]=-2;
644        if (colUpper[iColumn]>colLower[iColumn]) {
645          CoinBigIndex j;
646          double sum=0.0;
647          int iSOS=-1;
648          int n=0;
649          for (j = columnStart[iColumn];
650               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
651            int iRow = row[j];
652            if (contribution[iRow]>=0.0) {
653              iSOS=iRow;
654              n++;
655            } else {
656              sum += fabs(element[j]);
657            }
658          }
659          if (n>1)
660            COIN_DETAIL_PRINT(printf("Too many SOS entries (%d) for column %d\n",
661                                     n,iColumn));
662          if (sum) {
663            assert (iSOS>=0);
664            contribution[iSOS] += sum;
665            whichRow[iSOS]++;
666            whichSet[iColumn]=iSOS;
667          } else {
668            whichSet[iColumn]=iSOS+numberRows;
669          }
670        }
671      }
672      int * chunk = new int [numberRows];
673      for (int i=0;i<numberRows;i++) {
674        chunk[i]=-1;
675        if (whichRow[i]) {
676          contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]);
677        } else {
678          contribution[i] = COIN_DBL_MAX;
679        }
680        whichRow[i]=i;
681      }
682      newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
683      double * saveLower = CoinCopyOfArray(colLower,numberColumns);
684      double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
685      CoinSort_2(contribution,contribution+numberRows,whichRow);
686      // Set do nothing solution
687      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
688        if(whichSet[iColumn]>=numberRows)
689          newSolver->setColLower(iColumn,1.0);
690      }
691      newSolver->resolve();
692      int nChunk = (nSOS+9)/10;
693      int nPass=0;
694      int inChunk=0;
695      for (int i=0;i<nSOS;i++) {
696        chunk[whichRow[i]]=nPass;
697        inChunk++;
698        if (inChunk==nChunk) {
699          inChunk=0;
700          // last two together
701          if (i+nChunk<nSOS)
702            nPass++;
703        }
704      }
705      // adjust
706      nPass++;
707      for (int iPass=0;iPass<nPass;iPass++) {
708        // fix last chunk and unfix this chunk
709        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
710          int iSOS = whichSet[iColumn];
711          if (iSOS>=0) {
712            if (iSOS>=numberRows)
713              iSOS-=numberRows;
714            if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) {
715              newSolver->setColLower(iColumn,1.0);
716            } else if (chunk[iSOS]==iPass) {
717              newSolver->setColLower(iColumn,saveLower[iColumn]);
718              newSolver->setColUpper(iColumn,saveUpper[iColumn]);
719            }
720          }
721        }
722        // solve
723        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
724                                         model_->getCutoff(), "CbcHeuristicRENS");
725        if (returnCode < 0) {
726            returnCode = 0; // returned on size
727            break;
728        } else if ((returnCode&1)==0) {
729          // no good
730          break;
731        }
732      }
733      if ((returnCode&2) != 0) {
734        // could add cut
735        returnCode &= ~2;
736      }
737      delete [] chunk;
738      delete [] saveLower;
739      delete [] saveUpper;
740      delete [] whichRow;
741      delete [] whichSet;
742      delete [] contribution;
743      delete newSolver;
744      return returnCode;
745    }
746   
747    double primalTolerance;
748    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
749
750    int i;
751    int numberTightened = 0;
752    int numberAtBound = 0;
753    int numberContinuous = numberColumns - numberIntegers;
754    /*
755      0 - allow fixing
756      1 - don't allow fixing
757     */
758    char * marked = new char [numberColumns];
759    memset(marked,0,numberColumns);
760    if (bestSolution) {
761      for (i = 0; i < numberIntegers; i++) {
762        int iColumn = integerVariable[i];
763        double value = currentSolution[iColumn];
764        double lower = colLower[iColumn];
765        double upper = colUpper[iColumn];
766        value = CoinMax(value, lower);
767        value = CoinMin(value, upper);
768        if (fabs(bestSolution[iColumn]-value)>0.999) 
769          marked[iColumn]=1;
770      }
771    }
772    if ((rensType_&(64|128))!=0&&model_->objects()) {
773      int lowPriority=-COIN_INT_MAX;
774      int highPriority=COIN_INT_MAX;
775      for (i = 0; i < numberIntegers; i++) {
776        int priority=model_->priority(i);
777        lowPriority=CoinMax(lowPriority,priority);
778        highPriority=CoinMin(highPriority,priority);
779      }
780      if (highPriority<lowPriority) {
781        int keepPriority=((rensType_&64)!=0) ? highPriority : lowPriority; 
782        for (i = 0; i < numberIntegers; i++) {
783          int iColumn = integerVariable[i];
784          int priority=model_->priority(i);
785          if (priority==keepPriority)
786            marked[iColumn]=1;
787        }
788      }
789    }
790    for (i = 0; i < numberIntegers; i++) {
791        int iColumn = integerVariable[i];
792        double value = currentSolution[iColumn];
793        double lower = colLower[iColumn];
794        double upper = colUpper[iColumn];
795        value = CoinMax(value, lower);
796        value = CoinMin(value, upper);
797        double djValue=dj[iColumn]*direction;
798        bool dontFix=marked[iColumn]!=0;
799#define RENS_FIX_ONLY_LOWER
800#ifndef RENS_FIX_ONLY_LOWER
801        if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
802            value = floor(value + 0.5);
803            if (dontFix) {
804              continue;
805            }
806            if (value == lower || value == upper)
807                numberAtBound++;
808            newSolver->setColLower(iColumn, value);
809            newSolver->setColUpper(iColumn, value);
810            numberFixed++;
811        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
812            numberTightened++;
813            newSolver->setColLower(iColumn, floor(value));
814            newSolver->setColUpper(iColumn, ceil(value));
815        }
816#else
817        if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
818                floor(value + 0.5) == lower &&
819            djValue > djTolerance ) {
820            if (dontFix) {
821              continue;
822            }
823          value = floor(value + 0.5);
824          numberAtBound++;
825          newSolver->setColLower(iColumn, value);
826          newSolver->setColUpper(iColumn, value);
827          numberFixed++;
828        } else if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
829                floor(value + 0.5) == upper &&
830                   -djValue > djTolerance && (djTolerance > 0.0||type==2)) {
831          value = floor(value + 0.5);
832            if (dontFix) {
833              continue;
834            }
835          numberAtBound++;
836          newSolver->setColLower(iColumn, value);
837          newSolver->setColUpper(iColumn, value);
838          numberFixed++;
839        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0 &&
840                   djTolerance <0.0) {
841            numberTightened++;
842            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
843                value = floor(value + 0.5);
844                if (value < upper) {
845                    newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
846                    newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
847                } else {
848                    newSolver->setColLower(iColumn, upper - 1.0);
849                }
850            } else {
851                newSolver->setColLower(iColumn, floor(value));
852                newSolver->setColUpper(iColumn, ceil(value));
853            }
854        }
855#endif
856    }
857    delete [] marked;
858    delete [] dj;
859    if (numberFixed > numberIntegers / 5) {
860        if ( numberFixed < numberColumns / 5) {
861#define RENS_FIX_CONTINUOUS
862#ifdef RENS_FIX_CONTINUOUS
863            const double * colLower = newSolver->getColLower();
864            //const double * colUpper = newSolver->getColUpper();
865            int nAtLb = 0;
866            double sumDj = 0.0;
867            const double * dj = newSolver->getReducedCost();
868            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
869                if (!newSolver->isInteger(iColumn)) {
870                    double value = currentSolution[iColumn];
871                    if (value < colLower[iColumn] + 1.0e-8) {
872                        double djValue = dj[iColumn] * direction;
873                        nAtLb++;
874                        sumDj += djValue;
875                    }
876                }
877            }
878            if (nAtLb) {
879                // fix some continuous
880                double * sort = new double[nAtLb];
881                int * which = new int [nAtLb];
882                double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
883                int nFix2 = 0;
884                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
885                    if (!newSolver->isInteger(iColumn)) {
886                        double value = currentSolution[iColumn];
887                        if (value < colLower[iColumn] + 1.0e-8) {
888                            double djValue = dj[iColumn] * direction;
889                            if (djValue > threshold) {
890                                sort[nFix2] = -djValue;
891                                which[nFix2++] = iColumn;
892                            }
893                        }
894                    }
895                }
896                CoinSort_2(sort, sort + nFix2, which);
897                nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
898                for (int i = 0; i < nFix2; i++) {
899                    int iColumn = which[i];
900                    newSolver->setColUpper(iColumn, colLower[iColumn]);
901                }
902                delete [] sort;
903                delete [] which;
904#ifdef CLP_INVESTIGATE2
905                printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
906                       numberFixed, numberTightened, numberAtBound, nFix2);
907#endif
908            }
909#endif
910        }
911#ifdef COIN_DEVELOP
912        printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
913#endif
914        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
915                                         model_->getCutoff(), "CbcHeuristicRENS");
916        if (returnCode < 0 || returnCode == 0) {
917#ifdef RENS_FIX_CONTINUOUS
918            if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
919                const double * colLower = newSolver->getColLower();
920                //const double * colUpper = newSolver->getColUpper();
921                int nAtLb = 0;
922                double sumDj = 0.0;
923                const double * dj = newSolver->getReducedCost();
924                double direction = newSolver->getObjSense();
925                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
926                    if (!newSolver->isInteger(iColumn)) {
927                        double value = currentSolution[iColumn];
928                        if (value < colLower[iColumn] + 1.0e-8) {
929                            double djValue = dj[iColumn] * direction;
930                            nAtLb++;
931                            sumDj += djValue;
932                        }
933                    }
934                }
935                if (nAtLb) {
936                    // fix some continuous
937                    double * sort = new double[nAtLb];
938                    int * which = new int [nAtLb];
939                    double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
940                    int nFix2 = 0;
941                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
942                        if (!newSolver->isInteger(iColumn)) {
943                            double value = currentSolution[iColumn];
944                            if (value < colLower[iColumn] + 1.0e-8) {
945                                double djValue = dj[iColumn] * direction;
946                                if (djValue > threshold) {
947                                    sort[nFix2] = -djValue;
948                                    which[nFix2++] = iColumn;
949                                }
950                            }
951                        }
952                    }
953                    CoinSort_2(sort, sort + nFix2, which);
954                    nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
955                    for (int i = 0; i < nFix2; i++) {
956                        int iColumn = which[i];
957                        newSolver->setColUpper(iColumn, colLower[iColumn]);
958                    }
959                    delete [] sort;
960                    delete [] which;
961#ifdef CLP_INVESTIGATE2
962                    printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
963                           numberFixed, numberTightened, numberAtBound, nFix2);
964#endif
965                }
966                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
967                                                 model_->getCutoff(), "CbcHeuristicRENS");
968            }
969#endif
970            if (returnCode < 0 || returnCode == 0) {
971                // Do passes fixing up those >0.9 and
972                // down those < 0.05
973#define RENS_PASS 3
974              //#define KEEP_GOING
975#ifdef KEEP_GOING
976                double * saveLower = CoinCopyOfArray(colLower,numberColumns);
977                double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
978                bool badPass=false;
979                int nSolved=0;
980#endif
981                for (int iPass=0;iPass<RENS_PASS;iPass++) {
982                  int nFixed=0;
983                  int nFixedAlready=0;
984                  int nFixedContinuous=0;
985                  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
986                    if (colUpper[iColumn]>colLower[iColumn]) {
987                      if (newSolver->isInteger(iColumn)) {
988                        double value = currentSolution[iColumn];
989                        double fixTo = floor(value+0.1);
990                        if (fixTo>value || value-fixTo < 0.05) {
991                          // above 0.9 or below 0.05
992                          nFixed++;
993                          newSolver->setColLower(iColumn, fixTo);
994                          newSolver->setColUpper(iColumn, fixTo);
995                        }
996                      }
997                    } else if (newSolver->isInteger(iColumn)) {
998                      nFixedAlready++;
999                    } else {
1000                      nFixedContinuous++;
1001                    }
1002                  }
1003#ifdef CLP_INVESTIGATE2
1004                  printf("%d more integers fixed (total %d) plus %d continuous\n",
1005                         nFixed,nFixed+nFixedAlready,nFixedContinuous);
1006#endif
1007#ifdef KEEP_GOING
1008                  if (nFixed) {
1009                    newSolver->resolve();
1010                    if (!newSolver->isProvenOptimal()) {
1011                      badPass=true;
1012                      break;
1013                    } else {
1014                      nSolved++;
1015                      memcpy(saveLower,colLower,numberColumns*sizeof(double));
1016                      memcpy(saveUpper,colUpper,numberColumns*sizeof(double));
1017                    }
1018                  } else {
1019                    break;
1020                  }
1021#else
1022                  if (nFixed) {
1023                    newSolver->resolve();
1024                    if (!newSolver->isProvenOptimal()) {
1025                      returnCode=0;
1026                      break;
1027                    }
1028                    returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
1029                                                     model_->getCutoff(), "CbcHeuristicRENS");
1030                  } else {
1031                    returnCode=0;
1032                  }
1033                  if (returnCode>=0)
1034                    break;
1035                }
1036                if (returnCode < 0) 
1037                returnCode = 0; // returned on size
1038#endif
1039            }
1040#ifdef KEEP_GOING
1041            if (badPass) {
1042              newSolver->setColLower(saveLower);
1043              newSolver->setColUpper(saveUpper);
1044              newSolver->resolve();
1045            }
1046            delete [] saveLower;
1047            delete [] saveUpper;
1048            if (nSolved)
1049              returnCode = 
1050                smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
1051                                    model_->getCutoff(), "CbcHeuristicRENS");
1052            else
1053              returnCode=0;
1054        }
1055#endif
1056        }
1057        //printf("return code %d",returnCode);
1058        if ((returnCode&2) != 0) {
1059            // could add cut
1060            returnCode &= ~2;
1061#ifdef COIN_DEVELOP
1062            if (!numberTightened && numberFixed == numberAtBound)
1063                printf("could add cut with %d elements\n", numberFixed);
1064#endif
1065        } else {
1066            //printf("\n");
1067        }
1068    }
1069    //delete [] whichRow;
1070    //delete [] contribution;
1071    delete newSolver;
1072    fractionSmall_ = saveFractionSmall;
1073    return returnCode;
1074}
1075// update model
1076void CbcHeuristicRENS::setModel(CbcModel * model)
1077{
1078    model_ = model;
1079}
1080
Note: See TracBrowser for help on using the repository browser.