Changeset 1585


Ignore:
Timestamp:
Jan 11, 2011 2:04:34 PM (8 years ago)
Author:
forrest
Message:

add some more heuristics

Location:
trunk/Cbc/src
Files:
7 edited

Legend:

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

    r1582 r1585  
    838838    if (solver->isProvenOptimal()) {
    839839        CglPreProcess process;
     840        if ((model_->moreSpecialOptions()&65536)!=0)
     841          process.setOptions(2+4+8); // no cuts
    840842        /* Do not try and produce equality cliques and
    841843           do up to 2 passes (normally) 5 if restart */
  • trunk/Cbc/src/CbcHeuristicGreedy.cpp

    r1573 r1585  
    10641064    }
    10651065    double offset2 = 0.0;
     1066    char * sos = new char [numberRows];
    10661067    for (int iRow = 0;iRow < numberRows; iRow++) {
     1068      sos[iRow]=0;
     1069      if (rhs[iRow]<0.0) {
     1070        sos[iRow]=1;
     1071        rhs[iRow]=1.0;
     1072      }
    10671073      if( slackCost[iRow] == 1.0e30) {
    10681074        slackCost[iRow]=0.0;
    10691075      } else {
    10701076        offset2 += slackCost[iRow];
    1071         rhs[iRow] = -2.0;
     1077        sos[iRow] = 2;
    10721078      }
    10731079    }
     
    10951101    double * rowActivity = new double[numberRows];
    10961102    memset(rowActivity, 0, numberRows*sizeof(double));
    1097     if ((algorithm_&2)==0) {
     1103    if ((algorithm_&(2|4))==0) {
    10981104      // get solution as small as possible
    10991105      for (iColumn = 0; iColumn < numberColumns; iColumn++)
     
    11151121      }
    11161122    }
     1123    // get row activity
     1124    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1125      CoinBigIndex j;
     1126      double value = newSolution[iColumn];
     1127      for (j = columnStart[iColumn];
     1128           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1129        int iRow = row[j];
     1130        rowActivity[iRow] += value * element[j];
     1131      }
     1132    }
    11171133    double * contribution = new double [numberColumns];
    11181134    int * which = new int [numberColumns];
     
    11231139      double forSort = 0.0;
    11241140      bool hasSlack=false;
     1141      bool willFit=true;
    11251142      newSolutionValue += value * cost;
    11261143      for (j = columnStart[iColumn];
    11271144           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    11281145        int iRow = row[j];
    1129         rowActivity[iRow] += value * element[j];
    1130         if (rhs[iRow] >0.0) {
    1131           forSort += element[j];
    1132         } else  if (rhs[iRow] == -2.0) {
     1146        if (sos[iRow] == 2)
    11331147          hasSlack = true;
    1134         }
     1148        forSort += element[j];
     1149        double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
     1150        if (gap<element[j])
     1151          willFit = false;
    11351152      }
    1136       if (forSort && value == 0.0 && columnUpper[iColumn]) {
     1153      bool isSlack = hasSlack && (columnLength[iColumn]==1);
     1154      if ((algorithm_&4)!=0)
     1155        forSort=1.0;
     1156      // Use smallest cost if will fit
     1157      if (willFit && hasSlack &&
     1158          value == 0.0 && columnUpper[iColumn]) {
    11371159        if (hasSlack) {
    1138           if (cost>=0.0) {
     1160          if (cost>0.0) {
    11391161            forSort = 2.0e30;
     1162          } else if (cost==0.0) {
     1163            if (!isSlack)
     1164              forSort = 1.0e29;
     1165            else
     1166              forSort = 1.0e28;
    11401167          } else {
    11411168            forSort = cost/forSort;
     
    11631190           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    11641191        int iRow = row[j];
    1165         if (rhs[iRow]<0.0&&rowActivity[iRow]) {
     1192        if (sos[iRow]&&rowActivity[iRow]) {
    11661193          possible = false;
    1167         } else if (rhs[iRow]>=0.0) {
     1194        } else {
    11681195          double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
    11691196          if (gap<element[j])
     
    11831210      }
    11841211    }
     1212    delete [] sos;
    11851213    if (newSolutionValue < solutionValue) {
    11861214        // check feasible
  • trunk/Cbc/src/CbcHeuristicGreedy.hpp

    r1573 r1585  
    239239       1 bit - use current model, otherwise original
    240240       2 - use current solution as starting point, otherwise pure greedy
    241        4 - use duals to modify greedy
    242        8 - use duals on GUB/SOS in special way
     241       4 - as 2 but use merit not merit/size
     242       8 - use duals to modify greedy
     243       16 - use duals on GUB/SOS in special way
    243244    */
    244245    inline int algorithm() const {
  • trunk/Cbc/src/CbcHeuristicRENS.cpp

    r1582 r1585  
    9999    const double * currentSolution = newSolver->getColSolution();
    100100    int type = rensType_&15;
    101     if (type<11)
     101    if (type<12)
    102102      newSolver->resolve();
    103103    double direction = newSolver->getObjSense();
     
    108108    double tolerance;
    109109    newSolver->getDblParam(OsiDualTolerance, tolerance) ;
    110     if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<11) {
     110    if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<12) {
    111111      gap += 100.0 * tolerance;
    112112      int nFix = newSolver->reducedCostFix(gap);
     
    118118          << CoinMessageEol;
    119119      }
    120     } else if (type<11) {
     120    } else if (type<12) {
    121121      return 0; // finished?
    122122    }
     
    146146        delete basis;
    147147      }
    148     } else if (type>=5&&type<=11) {
     148    } else if (type>=5&&type<=12) {
    149149      /* 5 fix sets at one
    150150         6 fix on dj but leave unfixed SOS slacks
     
    153153         9 as 8 but only fix all at zero if just one in set nonzero
    154154         10 fix all "stable" ones
    155          11 layered approach
     155         11 fix all "stable" ones - approach 2
     156         12 layered approach
    156157      */
    157158      // SOS type fixing
    158       bool fixSets = (type==5)||(type==7)||(type==10);
     159      bool fixSets = (type==5)||(type==7)||(type==10)||(type==11);
    159160      CoinWarmStartBasis * basis =
    160161        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     
    250251        double * sort = new double [numberRows];
    251252        const double * pi = newSolver->getRowPrice();
    252         if (type==11) {
     253        if (type==12) {
    253254          contribution = new double [numberRows];
    254255          for (int i=0;i<numberRows;i++) {
     
    407408            numberFixed=numberColumns;
    408409            djTolerance = 1.0e30;
     410          } else if (type==11) {
     411            double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
     412            char * mark = new char [numberColumns];
     413            char * nonzero = new char [numberColumns];
     414            // save basis and solution
     415            CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(newSolver->getWarmStart()) ;
     416            assert(basis != NULL);
     417            double * saveSolution =
     418              CoinCopyOfArray(newSolver->getColSolution(),
     419                              numberColumns);
     420            double factors[] = {1.1,1.05,1.01,0.98};
     421            int nPass = (sizeof(factors)/sizeof(double))-1;
     422            double factor=factors[0];
     423            double proportion = fractionSmall_;
     424            fractionSmall_ = 0.5;
     425            // loosen up
     426            for (int i=0;i<numberRows;i++) {
     427              if (bestDj[i]>=1.0e30) {
     428                newSolver->setRowUpper(i,factor*saveUpper[i]);
     429              }
     430            }
     431            bool takeHint;
     432            OsiHintStrength strength;
     433            newSolver->getHintParam(OsiDoDualInResolve, takeHint, strength);
     434            newSolver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
     435            newSolver->resolve();
     436            newSolver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
     437            const double * solution = newSolver->getColSolution();
     438            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     439              mark[iColumn]=0;
     440              nonzero[iColumn]=0;
     441              if (colUpper[iColumn]>colLower[iColumn]&&
     442                  solution[iColumn]>0.9999)
     443                mark[iColumn]=1;
     444              else if (solution[iColumn]>0.00001)
     445                nonzero[iColumn]=1;
     446            }
     447            int nCheck=2;
     448            for (int iPass=0;iPass<nPass;iPass++) {
     449              // smaller
     450              factor = factors[iPass+1];
     451              for (int i=0;i<numberRows;i++) {
     452                if (bestDj[i]>=1.0e30) {
     453                  newSolver->setRowUpper(i,saveUpper[i]*factor);
     454                }
     455              }
     456              newSolver->resolve();
     457              if (newSolver->isProvenOptimal()) {
     458                nCheck++;
     459                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     460                  if (colUpper[iColumn]>colLower[iColumn]&&
     461                      solution[iColumn]>0.9999)
     462                    mark[iColumn]++;
     463                  else if (solution[iColumn]>0.00001)
     464                    nonzero[iColumn]++;
     465                }
     466              }
     467            }
     468            // correct values
     469            for (int i=0;i<numberRows;i++) {
     470              if (bestDj[i]>=1.0e30) {
     471                newSolver->setRowUpper(i,saveUpper[i]);
     472              }
     473            }
     474            newSolver->setColSolution(saveSolution);
     475            delete [] saveSolution;
     476            newSolver->setWarmStart(basis);
     477            delete basis ;
     478            newSolver->setHintParam(OsiDoDualInResolve, takeHint, strength);
     479            newSolver->resolve();
     480            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     481              if (colUpper[iColumn]>colLower[iColumn]&&
     482                  solution[iColumn]>0.9999)
     483                mark[iColumn]++;
     484              else if (solution[iColumn]>0.00001)
     485                nonzero[iColumn]++;
     486            }
     487            int nFixed=0;
     488            int numberSetsToFix = static_cast<int>(nSOS*(1.0-proportion));
     489            int * mixed = new int[numberRows];
     490            memset(mixed,0,numberRows*sizeof(int));
     491            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     492              if (colUpper[iColumn]>colLower[iColumn]) {
     493                int iSOS=-1;
     494                for (CoinBigIndex j = columnStart[iColumn];
     495                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     496                  int iRow = row[j];
     497                  if (bestDj[iRow]<1.0e25) {
     498                    iSOS=iRow;
     499                    break;
     500                  }
     501                }
     502                if (iSOS>=0) {
     503                  int numberTimesAtOne = mark[iColumn];
     504                  int numberTimesNonZero = nonzero[iColumn]+
     505                    numberTimesAtOne;
     506                  if (numberTimesAtOne<nCheck&&
     507                      numberTimesNonZero) {
     508                    mixed[iSOS]+=
     509                      CoinMin(numberTimesNonZero,
     510                              nCheck-numberTimesNonZero);
     511                  }
     512                }
     513              }
     514            }
     515            int nFix=0;
     516            for (int i=0;i<numberRows;i++) {
     517              if (bestDj[i]<1.0e25) {
     518                sort[nFix] = -bestDj[i]+1.0e8*mixed[i];
     519                mixed[nFix++]=i;
     520              }
     521            }
     522            CoinSort_2(sort,sort+nFix,mixed);
     523            nFix = CoinMin(nFix,numberSetsToFix);
     524            memset(sort,0,sizeof(double)*numberRows);
     525            for (int i=0;i<nFix;i++)
     526              sort[mixed[i]]=1.0;
     527            delete [] mixed;
     528            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     529              if (colUpper[iColumn]>colLower[iColumn]) {
     530                if (solution[iColumn]>0.9999) {
     531                  int iSOS=-1;
     532                  for (CoinBigIndex j = columnStart[iColumn];
     533                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     534                    int iRow = row[j];
     535                    if (bestDj[iRow]<1.0e25) {
     536                      iSOS=iRow;
     537                      break;
     538                    }
     539                  }
     540                  if (iSOS>=0&&sort[iSOS]) {
     541                    newSolver->setColLower(iColumn,1.0);
     542                    nFixed++;
     543                  }
     544                }
     545              }
     546            }
     547            char line[100];
     548            sprintf(line,"Heuristic %s fixed %d to one",
     549                    heuristicName(),
     550                    nFixed);
     551            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     552              << line
     553              << CoinMessageEol;
     554            delete [] mark;
     555            delete [] nonzero;
     556            delete [] saveUpper;
     557            numberFixed=numberColumns;
     558            djTolerance = 1.0e30;
    409559          }
    410560        }
     
    433583      djTolerance = CoinMax(sort[last],1.0e-5);
    434584      delete [] sort;
    435     } else if (type==11) {
     585    } else if (type==12) {
    436586      // Do layered in a different way
    437587      int numberRows = solver->getNumRows();
  • trunk/Cbc/src/CbcModel.cpp

    r1582 r1585  
    1038510385        if (continuousSolver_)
    1038610386            solver_ = continuousSolver_;
     10387        // save basis and solution
     10388        CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
     10389        assert(basis != NULL);
     10390        double * saveSolution = CoinCopyOfArray(solver_->getColSolution(),
     10391                                                solver_->getNumCols());
    1038710392        // move solution to continuous copy
    1038810393        solver_->setColSolution(solution);
     
    1069210697        delete [] saveUpper;
    1069310698
     10699        solver_->setColSolution(saveSolution);
     10700        delete [] saveSolution;
     10701        solver_->setWarmStart(basis);
     10702        delete basis ;
    1069410703        /*
    1069510704          Restore the usual solver.
  • trunk/Cbc/src/CbcModel.hpp

    r1582 r1585  
    17361736        16384 gomory
    17371737        32768 more heuristics in sub trees
     1738        65536 no cuts in preprocessing
    17381739    */
    17391740    inline void setMoreSpecialOptions(int value) {
  • trunk/Cbc/src/CbcSolver.cpp

    r1573 r1585  
    35883588                                        process.setOptions(1);
    35893589                                    // Add in generators
    3590                                     process.addCutGenerator(&generator1);
     3590                                    if ((model_.moreSpecialOptions()&65536)==0)
     3591                                      process.addCutGenerator(&generator1);
    35913592                                    int translate[] = {9999, 0, 0, -3, 2, 3, -2, 9999, 4, 5};
    35923593                                    process.passInMessageHandler(babModel_->messageHandler());
     
    36913692                                            osiclp->getModelPtr()->setPerturbation(52); // try less
    36923693#endif
     3694                                        if ((model_.moreSpecialOptions()&65536)!=0)
     3695                                          process.setOptions(2+4+8); // no cuts
    36933696                                        solver2 = process.preProcessNonDefault(*saveSolver, translate[preProcess], numberPasses,
    36943697                                                                               tunePreProcess);
     
    75057508                            // get next field
    75067509                            field = CoinReadGetString(argc, argv);
     7510                            bool append = false;
     7511                            if (field == "append$") {
     7512                              field = "$";
     7513                              append = true;
     7514                            }
    75077515                            if (field == "$") {
    75087516                                field = parameters_[iParam].stringValue();
     
    75477555                                    fileName = directory + field;
    75487556                                }
    7549                                 fp = fopen(fileName.c_str(), "w");
     7557                                if (!append)
     7558                                  fp = fopen(fileName.c_str(), "w");
     7559                                else
     7560                                  fp = fopen(fileName.c_str(), "a");
    75507561                            }
    75517562                            if (fp) {
    75527563#ifndef CBC_OTHER_SOLVER
    7553                                 if (printMode != 5) {
     7564                                if (printMode < 5) {
    75547565                                    // Write solution header (suggested by Luigi Poderico)
    75557566                                    lpSolver->computeObjectiveValue(false);
     
    77067717                                    delete [] newMasks;
    77077718                                }
     7719                                if (printMode > 5) {
     7720                                  ClpSimplex * solver = clpSolver->getModelPtr();
     7721                                  int numberColumns = solver->numberColumns();
     7722                                  // column length unless rhs ranging
     7723                                  int number = numberColumns;
     7724                                  switch (printMode) {
     7725                                    // bound ranging
     7726                                  case 6:
     7727                                    fprintf(fp,"Bound ranging");
     7728                                    break;
     7729                                    // rhs ranging
     7730                                  case 7:
     7731                                    fprintf(fp,"Rhs ranging");
     7732                                    number = numberRows;
     7733                                    break;
     7734                                    // objective ranging
     7735                                  case 8:
     7736                                    fprintf(fp,"Objective ranging");
     7737                                    break;
     7738                                  }
     7739                                  if (lengthName)
     7740                                    fprintf(fp,",name");
     7741                                  fprintf(fp,",increase,variable,decrease,variable\n");
     7742                                  int * which = new int [ number];
     7743                                  if (printMode != 7) {
     7744                                    if (!doMask) {
     7745                                      for (int i = 0; i < number;i ++)
     7746                                        which[i]=i;
     7747                                    } else {
     7748                                      int n = 0;
     7749                                      for (int i = 0; i < number;i ++) {
     7750                                        if (maskMatches(maskStarts,masks,columnNames[i]))
     7751                                          which[n++]=i;
     7752                                      }
     7753                                      if (n) {
     7754                                        number=n;
     7755                                      } else {
     7756                                        printf("No names match - doing all\n");
     7757                                        for (int i = 0; i < number;i ++)
     7758                                          which[i]=i;
     7759                                      }
     7760                                    }
     7761                                  } else {
     7762                                    if (!doMask) {
     7763                                    for (int i = 0; i < number;i ++)
     7764                                      which[i]=i+numberColumns;
     7765                                    } else {
     7766                                      int n = 0;
     7767                                      for (int i = 0; i < number;i ++) {
     7768                                        if (maskMatches(maskStarts,masks,rowNames[i]))
     7769                                          which[n++]=i+numberColumns;
     7770                                      }
     7771                                      if (n) {
     7772                                        number=n;
     7773                                      } else {
     7774                                        printf("No names match - doing all\n");
     7775                                        for (int i = 0; i < number;i ++)
     7776                                          which[i]=i+numberColumns;
     7777                                      }
     7778                                    }
     7779                                  }
     7780                                  double * valueIncrease = new double [ number];
     7781                                  int * sequenceIncrease = new int [ number];
     7782                                  double * valueDecrease = new double [ number];
     7783                                  int * sequenceDecrease = new int [ number];
     7784                                  switch (printMode) {
     7785                                    // bound or rhs ranging
     7786                                  case 6:
     7787                                  case 7:
     7788                                    solver->primalRanging(numberRows,
     7789                                                                 which, valueIncrease, sequenceIncrease,
     7790                                                                 valueDecrease, sequenceDecrease);
     7791                                    break;
     7792                                    // objective ranging
     7793                                  case 8:
     7794                                    solver->dualRanging(number,
     7795                                                               which, valueIncrease, sequenceIncrease,
     7796                                                               valueDecrease, sequenceDecrease);
     7797                                    break;
     7798                                  }
     7799                                  for (int i = 0; i < number; i++) {
     7800                                    int iWhich = which[i];
     7801                                    fprintf(fp, "%d,", (iWhich<numberColumns) ? iWhich : iWhich-numberColumns);
     7802                                    if (lengthName) {
     7803                                      const char * name = (printMode==7) ? rowNames[iWhich-numberColumns].c_str() : columnNames[iWhich].c_str();
     7804                                      fprintf(fp,"%s,",name);
     7805                                    }
     7806                                    if (valueIncrease[i]<1.0e30) {
     7807                                      fprintf(fp, "%.10g,", valueIncrease[i]);
     7808                                      int outSequence = sequenceIncrease[i];
     7809                                      if (outSequence<numberColumns) {
     7810                                        if (lengthName)
     7811                                          fprintf(fp,"%s,",columnNames[outSequence].c_str());
     7812                                        else
     7813                                          fprintf(fp,"C%7.7d,",outSequence);
     7814                                      } else {
     7815                                        outSequence -= numberColumns;
     7816                                        if (lengthName)
     7817                                          fprintf(fp,"%s,",rowNames[outSequence].c_str());
     7818                                        else
     7819                                          fprintf(fp,"R%7.7d,",outSequence);
     7820                                      }
     7821                                    } else {
     7822                                      fprintf(fp,"1.0e100,,");
     7823                                    }
     7824                                    if (valueDecrease[i]<1.0e30) {
     7825                                      fprintf(fp, "%.10g,", valueDecrease[i]);
     7826                                      int outSequence = sequenceDecrease[i];
     7827                                      if (outSequence<numberColumns) {
     7828                                        if (lengthName)
     7829                                          fprintf(fp,"%s",columnNames[outSequence].c_str());
     7830                                        else
     7831                                          fprintf(fp,"C%7.7d",outSequence);
     7832                                      } else {
     7833                                        outSequence -= numberColumns;
     7834                                        if (lengthName)
     7835                                          fprintf(fp,"%s",rowNames[outSequence].c_str());
     7836                                        else
     7837                                          fprintf(fp,"R%7.7d",outSequence);
     7838                                      }
     7839                                    } else {
     7840                                      fprintf(fp,"1.0e100,");
     7841                                    }
     7842                                    fprintf(fp,"\n");
     7843                                  }
     7844                                  if (fp != stdout)
     7845                                    fclose(fp);
     7846                                  delete [] which;
     7847                                  delete [] valueIncrease;
     7848                                  delete [] sequenceIncrease;
     7849                                  delete [] valueDecrease;
     7850                                  delete [] sequenceDecrease;
     7851                                  if (masks) {
     7852                                    delete [] maskStarts;
     7853                                    for (int i = 0; i < maxMasks; i++)
     7854                                      delete [] masks[i];
     7855                                    delete [] masks;
     7856                                  }
     7857                                  break;
     7858                                }
    77087859                                if (printMode > 2 && printMode < 5) {
    77097860                                    for (iRow = 0; iRow < numberRows; iRow++) {
     
    87368887  Source code changes so up to 2.0
    87378888*/
    8738 
    8739 
Note: See TracChangeset for help on using the changeset viewer.