Changeset 1582 for trunk


Ignore:
Timestamp:
Jan 7, 2011 12:16:00 PM (9 years ago)
Author:
forrest
Message:

add some sos heuristics

Location:
trunk/Cbc/src
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcHeuristic.cpp

    r1573 r1582  
    703703    int returnCode = 1;
    704704    int saveModelOptions = model_->specialOptions();
    705     assert ((saveModelOptions&2048) == 0);
     705    //assert ((saveModelOptions&2048) == 0);
    706706    model_->setSpecialOptions(saveModelOptions | 2048);
    707707    {
     
    924924                    model.setMaximumNodes(numberNodes);
    925925                    model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
     926                    if ((saveModelOptions&2048) == 0)
     927                      model.setMoreSpecialOptions(model_->moreSpecialOptions());
    926928                    // Lightweight
    927929                    CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
  • trunk/Cbc/src/CbcHeuristicRENS.cpp

    r1573 r1582  
    9090        return 0;
    9191    numberTries_++;
     92    double saveFractionSmall=fractionSmall_;
    9293    OsiSolverInterface * solver = model_->solver();
    9394
     
    9596    const int * integerVariable = model_->integerVariable();
    9697
    97     const double * currentSolution = solver->getColSolution();
    9898    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
    99     newSolver->resolve();
     99    const double * currentSolution = newSolver->getColSolution();
     100    int type = rensType_&15;
     101    if (type<11)
     102      newSolver->resolve();
    100103    double direction = newSolver->getObjSense();
    101104    double cutoff=model_->getCutoff();
    102     //newSolver->getDblParam(OsiDualObjectiveLimit, cutoff);
     105    newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
    103106    //cutoff *= direction;
    104107    double gap = cutoff - newSolver->getObjValue() * direction ;
    105108    double tolerance;
    106109    newSolver->getDblParam(OsiDualTolerance, tolerance) ;
    107     if (gap > 0.0 || !newSolver->isProvenOptimal()) {
     110    if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<11) {
    108111      gap += 100.0 * tolerance;
    109112      int nFix = newSolver->reducedCostFix(gap);
     
    115118          << CoinMessageEol;
    116119      }
    117     } else {
     120    } else if (type<11) {
    118121      return 0; // finished?
    119122    }
    120123    int numberColumns = solver->getNumCols();
    121124    double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
    122     int type = rensType_&15;
    123125    double djTolerance = (type!=1) ? -1.0e30 : 1.0e-4;
    124126    const double * colLower = newSolver->getColLower();
    125127    const double * colUpper = newSolver->getColUpper();
     128    double * contribution = NULL;
    126129    int numberFixed = 0;
    127130    if (type==3) {
     
    143146        delete basis;
    144147      }
    145     } else if (type>=5&&type<=10) {
     148    } else if (type>=5&&type<=11) {
    146149      /* 5 fix sets at one
    147150         6 fix on dj but leave unfixed SOS slacks
     
    149152         8 fix all at zero but leave unfixed SOS slacks
    150153         9 as 8 but only fix all at zero if just one in set nonzero
    151          10 as 7 but pi other way
     154         10 fix all "stable" ones
     155         11 layered approach
    152156      */
    153157      // SOS type fixing
     
    203207          }
    204208        }
     209        // Just leave one slack in each set
     210        {
     211          const double * objective = newSolver->getObjCoefficients();
     212          int * best = new int [numberRows];
     213          double * cheapest = new double[numberRows];
     214          for (int i=0;i<numberRows;i++) {
     215            best[i]=-1;
     216            cheapest[i]=COIN_DBL_MAX;
     217          }
     218          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     219            if (colUpper[iColumn]>colLower[iColumn]) {
     220              if (columnLength[iColumn]==1) {
     221                CoinBigIndex j = columnStart[iColumn];
     222                int iRow = row[j];
     223                if (bestDj[iRow]<1.0e30) {
     224                  double obj = direction*objective[iColumn];
     225                  if (obj<cheapest[iRow]) {
     226                    cheapest[iRow]=obj;
     227                    best[iRow]=iColumn;
     228                  }
     229                }
     230              }
     231            }
     232          }
     233          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     234            if (colUpper[iColumn]>colLower[iColumn]) {
     235              if (columnLength[iColumn]==1) {
     236                CoinBigIndex j = columnStart[iColumn];
     237                int iRow = row[j];
     238                if (bestDj[iRow]<1.0e30) {
     239                  if (best[iRow]!=-1&&iColumn!=best[iRow]) {
     240                    newSolver->setColUpper(iColumn,0.0);
     241                  }
     242                }
     243              }
     244            }
     245          }
     246          delete [] best;
     247          delete [] cheapest;
     248        }
    205249        int nSOS=0;
    206250        double * sort = new double [numberRows];
    207251        const double * pi = newSolver->getRowPrice();
     252        if (type==11) {
     253          contribution = new double [numberRows];
     254          for (int i=0;i<numberRows;i++) {
     255            if (bestDj[i]<1.0e30)
     256              contribution[i]=0.0;
     257            else
     258              contribution[i]=-1.0;
     259          }
     260        }
    208261        for (int i=0;i<numberRows;i++) {
    209262          if (bestDj[i]<1.0e30) {
     
    217270        }
    218271        if (10*nSOS>8*numberRows) {
    219           std::sort(sort,sort+nSOS);
    220           int last = static_cast<int>(nSOS*0.9*fractionSmall_);
    221           double tolerance = sort[last];
    222           for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    223             if (colUpper[iColumn]>colLower[iColumn]) {
    224               CoinBigIndex j;
    225               if (currentSolution[iColumn]<=1.0e-6||
    226                   currentSolution[iColumn]>=0.999999) {
    227                 if (fixSets) {
    228                   for (j = columnStart[iColumn];
    229                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    230                     int iRow = row[j];
    231                     double useDj;
    232                     if (type==5)
    233                       useDj = bestDj[iRow];
    234                     else if (type==7)
    235                       useDj= -fabs(pi[iRow]);
    236                     else
    237                       useDj= fabs(pi[iRow]);
    238                     if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
    239                       numberFixed++;
    240                       if (currentSolution[iColumn]<=1.0e-6)
    241                         newSolver->setColUpper(iColumn,0.0);
    242                       else if (currentSolution[iColumn]>=0.999999)
    243                         newSolver->setColLower(iColumn,1.0);
     272          if (type<10) {
     273            std::sort(sort,sort+nSOS);
     274            int last = static_cast<int>(nSOS*0.9*fractionSmall_);
     275            double tolerance = sort[last];
     276            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     277              if (colUpper[iColumn]>colLower[iColumn]) {
     278                CoinBigIndex j;
     279                if (currentSolution[iColumn]<=1.0e-6||
     280                    currentSolution[iColumn]>=0.999999) {
     281                  if (fixSets) {
     282                    for (j = columnStart[iColumn];
     283                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     284                      int iRow = row[j];
     285                      double useDj;
     286                      if (type==5)
     287                        useDj = bestDj[iRow];
     288                      else if (type==7)
     289                        useDj= -fabs(pi[iRow]);
     290                      else
     291                        useDj= fabs(pi[iRow]);
     292                      if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
     293                        numberFixed++;
     294                        if (currentSolution[iColumn]<=1.0e-6)
     295                          newSolver->setColUpper(iColumn,0.0);
     296                        else if (currentSolution[iColumn]>=0.999999)
     297                          newSolver->setColLower(iColumn,1.0);
     298                      }
     299                    }
     300                  } else if (columnLength[iColumn]==1) {
     301                    // leave more slacks
     302                    int iRow = row[columnStart[iColumn]];
     303                    if (bestDj[iRow]<1.0e30) {
     304                      // fake dj
     305                      dj[iColumn] *= 0.000001;
     306                    }
     307                  } else if (type==8||type==9) {
     308                    if (currentSolution[iColumn]<=1.0e-6) {
     309                      if (type==8) {
     310                        dj[iColumn] *= 1.0e6;
     311                      } else {
     312                        bool fix=false;
     313                        for (j = columnStart[iColumn];
     314                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     315                          int iRow = row[j];
     316                          if (bestDj[iRow]<1.0e25) {
     317                            fix=true;
     318                            break;
     319                          }
     320                        }
     321                        if (fix) {
     322                          dj[iColumn] *= 1.0e6;
     323                        }
     324                      }
     325                    } else {
     326                      dj[iColumn] *= 0.000001;
    244327                    }
    245328                  }
    246                 } else if (columnLength[iColumn]==1) {
    247                   // leave more slacks
    248                   int iRow = row[columnStart[iColumn]];
    249                   if (bestDj[iRow]<1.0e30) {
    250                     // fake dj
    251                     dj[iColumn] *= 0.000001;
    252                   }
    253                 } else if (type==8||type==9) {
    254                   if (currentSolution[iColumn]<=1.0e-6) {
    255                     if (type==8) {
    256                       dj[iColumn] *= 1.0e6;
    257                     } else {
    258                       bool fix=false;
    259                       for (j = columnStart[iColumn];
    260                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    261                         int iRow = row[j];
    262                         if (bestDj[iRow]<1.0e25) {
    263                           fix=true;
    264                           break;
    265                         }
    266                       }
    267                       if (fix) {
    268                         dj[iColumn] *= 1.0e6;
    269                       }
    270                     }
    271                   } else {
    272                     dj[iColumn] *= 0.000001;
    273                   }
    274329                }
    275330              }
    276331            }
    277           }
    278           if (fixSets)
     332            if (fixSets)
     333              djTolerance = 1.0e30;
     334          } else if (type==10) {
     335            double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
     336            char * mark = new char [numberColumns];
     337            char * nonzero = new char [numberColumns];
     338            double factor=CoinMax(1.000001,fractionSmall_);
     339            fractionSmall_ = 0.5;
     340            // loosen up
     341            for (int i=0;i<numberRows;i++) {
     342              if (bestDj[i]>=1.0e30) {
     343                newSolver->setRowUpper(i,factor*saveUpper[i]);
     344              }
     345            }
     346            newSolver->resolve();
     347            const double * solution = newSolver->getColSolution();
     348            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     349              mark[iColumn]=0;
     350              nonzero[iColumn]=0;
     351              if (colUpper[iColumn]>colLower[iColumn]&&
     352                  solution[iColumn]>0.9999)
     353                mark[iColumn]=1;
     354              else if (solution[iColumn]>0.00001)
     355                nonzero[iColumn]=1;
     356            }
     357            // slightly small
     358            for (int i=0;i<numberRows;i++) {
     359              if (bestDj[i]>=1.0e30) {
     360                newSolver->setRowUpper(i,saveUpper[i]*0.9999);
     361              }
     362            }
     363            newSolver->resolve();
     364            int nCheck=2;
     365            if (newSolver->isProvenOptimal()) {
     366              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     367                if (colUpper[iColumn]>colLower[iColumn]&&
     368                    solution[iColumn]>0.9999)
     369                  mark[iColumn]++;
     370                else if (solution[iColumn]>0.00001)
     371                  nonzero[iColumn]=1;
     372              }
     373            } else {
     374              nCheck=1;
     375            }
     376            // correct values
     377            for (int i=0;i<numberRows;i++) {
     378              if (bestDj[i]>=1.0e30) {
     379                newSolver->setRowUpper(i,saveUpper[i]);
     380              }
     381            }
     382            newSolver->resolve();
     383            int nFixed=0;
     384            int nFixedToZero=0;
     385            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     386              if (colUpper[iColumn]>colLower[iColumn]) {
     387                if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) {
     388                  newSolver->setColLower(iColumn,1.0);
     389                  nFixed++;
     390                } else if (!mark[iColumn]&&!nonzero[iColumn]&&
     391                           columnLength[iColumn]>1&&solution[iColumn]<0.00001) {
     392                  newSolver->setColUpper(iColumn,0.0);
     393                  nFixedToZero++;
     394                }
     395              }
     396            }
     397            char line[100];
     398            sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)",
     399                    heuristicName(),
     400                    nFixed,nFixedToZero);
     401            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     402              << line
     403              << CoinMessageEol;
     404            delete [] mark;
     405            delete []nonzero;
     406            delete [] saveUpper;
     407            numberFixed=numberColumns;
    279408            djTolerance = 1.0e30;
     409          }
    280410        }
    281411        delete basis;
    282412        delete [] sort;
    283413        delete [] bestDj;
     414        if (10*nSOS<=8*numberRows) {
     415          // give up
     416          delete [] contribution;
     417          delete newSolver;
     418          return 0;
     419        }
    284420      }
    285421    }
    286422    // Do dj to get right number
    287     if (type==4||type==6||type>7) {
     423    if (type==4||type==6||(type>7&&type<10)) {
    288424      double * sort = new double [numberColumns];
    289425      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     
    297433      djTolerance = CoinMax(sort[last],1.0e-5);
    298434      delete [] sort;
     435    } else if (type==11) {
     436      // Do layered in a different way
     437      int numberRows = solver->getNumRows();
     438      // Column copy
     439      const CoinPackedMatrix * matrix = newSolver->getMatrixByCol();
     440      const double * element = matrix->getElements();
     441      const int * row = matrix->getIndices();
     442      const CoinBigIndex * columnStart = matrix->getVectorStarts();
     443      const int * columnLength = matrix->getVectorLengths();
     444      int * whichRow = new int[numberRows];
     445      int * whichSet = new int [numberColumns];
     446      int nSOS=0;
     447      for (int i=0;i<numberRows;i++) {
     448        whichRow[i]=0;
     449        if (!contribution[i])
     450          nSOS++;
     451      }
     452      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     453        whichSet[iColumn]=-2;
     454        if (colUpper[iColumn]>colLower[iColumn]) {
     455          CoinBigIndex j;
     456          double sum=0.0;
     457          int iSOS=-1;
     458          int n=0;
     459          for (j = columnStart[iColumn];
     460               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     461            int iRow = row[j];
     462            if (contribution[iRow]>=0.0) {
     463              iSOS=iRow;
     464              n++;
     465            } else {
     466              sum += fabs(element[j]);
     467            }
     468          }
     469          if (n>1)
     470            printf("Too many SOS entries (%d) for column %d\n",
     471                   n,iColumn);
     472          if (sum) {
     473            assert (iSOS>=0);
     474            contribution[iSOS] += sum;
     475            whichRow[iSOS]++;
     476            whichSet[iColumn]=iSOS;
     477          } else {
     478            whichSet[iColumn]=iSOS+numberRows;
     479          }
     480        }
     481      }
     482      int * chunk = new int [numberRows];
     483      for (int i=0;i<numberRows;i++) {
     484        chunk[i]=-1;
     485        if (whichRow[i]) {
     486          contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]);
     487        } else {
     488          contribution[i] = COIN_DBL_MAX;
     489        }
     490        whichRow[i]=i;
     491      }
     492      newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
     493      double * saveLower = CoinCopyOfArray(colLower,numberColumns);
     494      double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
     495      CoinSort_2(contribution,contribution+numberRows,whichRow);
     496      // Set do nothing solution
     497      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     498        if(whichSet[iColumn]>=numberRows)
     499          newSolver->setColLower(iColumn,1.0);
     500      }
     501      newSolver->resolve();
     502      int nChunk = (nSOS+9)/10;
     503      int nPass=0;
     504      int inChunk=0;
     505      for (int i=0;i<nSOS;i++) {
     506        chunk[whichRow[i]]=nPass;
     507        inChunk++;
     508        if (inChunk==nChunk) {
     509          inChunk=0;
     510          // last two together
     511          if (i+nChunk<nSOS)
     512            nPass++;
     513        }
     514      }
     515      // adjust
     516      nPass++;
     517      for (int iPass=0;iPass<nPass;iPass++) {
     518        // fix last chunk and unfix this chunk
     519        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     520          int iSOS = whichSet[iColumn];
     521          if (iSOS>=0) {
     522            if (iSOS>=numberRows)
     523              iSOS-=numberRows;
     524            if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) {
     525              newSolver->setColLower(iColumn,1.0);
     526            } else if (chunk[iSOS]==iPass) {
     527              newSolver->setColLower(iColumn,saveLower[iColumn]);
     528              newSolver->setColUpper(iColumn,saveUpper[iColumn]);
     529            }
     530          }
     531        }
     532        // solve
     533        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
     534                                         model_->getCutoff(), "CbcHeuristicRENS");
     535        if (returnCode < 0) {
     536            returnCode = 0; // returned on size
     537            break;
     538        } else if ((returnCode&1)==0) {
     539          // no good
     540          break;
     541        }
     542      }
     543      if ((returnCode&2) != 0) {
     544        // could add cut
     545        returnCode &= ~2;
     546      }
     547      delete [] chunk;
     548      delete [] saveLower;
     549      delete [] saveUpper;
     550      delete [] whichRow;
     551      delete [] whichSet;
     552      delete [] contribution;
     553      delete newSolver;
     554      return returnCode;
    299555    }
    300556   
     
    490746        }
    491747    }
    492 
     748    //delete [] whichRow;
     749    //delete [] contribution;
    493750    delete newSolver;
     751    fractionSmall_ = saveFractionSmall;
    494752    return returnCode;
    495753}
  • trunk/Cbc/src/CbcModel.cpp

    r1573 r1582  
    1287112871                    if (handler_->logLevel()>1) {
    1287212872                      char line[100];
    12873                       sprintf(line,"Heuristic %s took %g seconds",
     12873                      sprintf(line,"Heuristic %s took %g seconds (%s)",
    1287412874                              heuristic_[i]->heuristicName(),
    12875                               getCurrentSeconds()-before);
     12875                              getCurrentSeconds()-before,
     12876                              ifSol ? "good" : "no good");
    1287612877                      handler_->message(CBC_GENERAL, messages_) <<
    1287712878                        line << CoinMessageEol ;
  • trunk/Cbc/src/CbcModel.hpp

    r1573 r1582  
    17351735        8192 slowly increase minimum drop
    17361736        16384 gomory
     1737        32768 more heuristics in sub trees
    17371738    */
    17381739    inline void setMoreSpecialOptions(int value) {
  • trunk/Cbc/src/CbcStrategy.cpp

    r1573 r1582  
    4242#include "CbcHeuristic.hpp"
    4343#include "CbcHeuristicLocal.hpp"
     44#include "CbcHeuristicRINS.hpp"
    4445
    4546// Default Constructor
     
    925926    if (!found)
    926927        model.addHeuristic(&heuristic1);
     928    if ((model.moreSpecialOptions()&32768)!=0) {
     929      // Allow join solutions
     930      CbcHeuristicLocal heuristic2(model);
     931      heuristic2.setHeuristicName("join solutions");
     932      //sheuristic2.setSearchType(1);
     933      found = false;
     934      for (iHeuristic = 0; iHeuristic < numberHeuristics; iHeuristic++) {
     935        CbcHeuristic * heuristic = model.heuristic(iHeuristic);
     936        CbcHeuristicLocal * cgl = dynamic_cast<CbcHeuristicLocal *>(heuristic);
     937        if (cgl) {
     938          found = true;
     939          break;
     940        }
     941      }
     942      if (!found)
     943        model.addHeuristic(&heuristic2);
     944      // Allow RINS
     945      CbcHeuristicRINS heuristic5(model);
     946      heuristic5.setHeuristicName("RINS");
     947      heuristic5.setFractionSmall(0.5);
     948      heuristic5.setDecayFactor(5.0);
     949      //heuristic5.setSearchType(1);
     950      found = false;
     951      for (iHeuristic = 0; iHeuristic < numberHeuristics; iHeuristic++) {
     952        CbcHeuristic * heuristic = model.heuristic(iHeuristic);
     953        CbcHeuristicLocal * cgl = dynamic_cast<CbcHeuristicLocal *>(heuristic);
     954        if (cgl) {
     955          found = true;
     956          break;
     957        }
     958      }
     959      if (!found)
     960        model.addHeuristic(&heuristic5);
     961    }
    927962}
    928963// Do printing stuff
Note: See TracChangeset for help on using the changeset viewer.