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

Last change on this file since 2104 was 2094, checked in by forrest, 5 years ago

for memory leaks and heuristics and some experimental stuff

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