Ignore:
Timestamp:
Dec 30, 2010 12:45:15 PM (9 years ago)
Author:
forrest
Message:

add some heuristic variants

File:
1 edited

Legend:

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

    r1286 r1564  
    861861}
    862862
    863 
     863// Default Constructor
     864CbcHeuristicGreedySOS::CbcHeuristicGreedySOS()
     865        : CbcHeuristic()
     866{
     867    originalRhs_ = NULL;
     868    // matrix  will automatically be empty
     869    originalNumberRows_ = 0;
     870    algorithm_ = 0;
     871    numberTimes_ = 100;
     872}
     873
     874// Constructor from model
     875CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(CbcModel & model)
     876        : CbcHeuristic(model)
     877{
     878    gutsOfConstructor(&model);
     879    algorithm_ = 2;
     880    numberTimes_ = 100;
     881    whereFrom_ = 1;
     882}
     883
     884// Destructor
     885CbcHeuristicGreedySOS::~CbcHeuristicGreedySOS ()
     886{
     887  delete [] originalRhs_;
     888}
     889
     890// Clone
     891CbcHeuristic *
     892CbcHeuristicGreedySOS::clone() const
     893{
     894    return new CbcHeuristicGreedySOS(*this);
     895}
     896// Guts of constructor from a CbcModel
     897void
     898CbcHeuristicGreedySOS::gutsOfConstructor(CbcModel * model)
     899{
     900    model_ = model;
     901    // Get a copy of original matrix
     902    assert(model->solver());
     903    if (model->solver()->getNumRows()) {
     904        matrix_ = *model->solver()->getMatrixByCol();
     905    }
     906    originalNumberRows_ = model->solver()->getNumRows();
     907    originalRhs_ = new double [originalNumberRows_];
     908}
     909// Create C++ lines to get to current state
     910void
     911CbcHeuristicGreedySOS::generateCpp( FILE * fp)
     912{
     913    CbcHeuristicGreedySOS other;
     914    fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n");
     915    fprintf(fp, "3  CbcHeuristicGreedySOS heuristicGreedySOS(*cbcModel);\n");
     916    CbcHeuristic::generateCpp(fp, "heuristicGreedySOS");
     917    if (algorithm_ != other.algorithm_)
     918        fprintf(fp, "3  heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_);
     919    else
     920        fprintf(fp, "4  heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_);
     921    if (numberTimes_ != other.numberTimes_)
     922        fprintf(fp, "3  heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_);
     923    else
     924        fprintf(fp, "4  heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_);
     925    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicGreedySOS);\n");
     926}
     927
     928// Copy constructor
     929CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(const CbcHeuristicGreedySOS & rhs)
     930        :
     931        CbcHeuristic(rhs),
     932        matrix_(rhs.matrix_),
     933        originalNumberRows_(rhs.originalNumberRows_),
     934        algorithm_(rhs.algorithm_),
     935        numberTimes_(rhs.numberTimes_)
     936{
     937  originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_);
     938}
     939
     940// Assignment operator
     941CbcHeuristicGreedySOS &
     942CbcHeuristicGreedySOS::operator=( const CbcHeuristicGreedySOS & rhs)
     943{
     944    if (this != &rhs) {
     945        CbcHeuristic::operator=(rhs);
     946        matrix_ = rhs.matrix_;
     947        originalNumberRows_ = rhs.originalNumberRows_;
     948        algorithm_ = rhs.algorithm_;
     949        numberTimes_ = rhs.numberTimes_;
     950        delete [] originalRhs_;
     951        originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_);
     952    }
     953    return *this;
     954}
     955// Returns 1 if solution, 0 if not
     956int
     957CbcHeuristicGreedySOS::solution(double & solutionValue,
     958                                  double * betterSolution)
     959{
     960    numCouldRun_++;
     961    if (!model_)
     962        return 0;
     963    // See if to do
     964    if (!when() || (when() == 1 && model_->phase() != 1))
     965        return 0; // switched off
     966    if (model_->getNodeCount() > numberTimes_)
     967        return 0;
     968    // See if at root node
     969    bool atRoot = model_->getNodeCount() == 0;
     970    int passNumber = model_->getCurrentPassNumber();
     971    if (atRoot && passNumber != 1)
     972        return 0;
     973    OsiSolverInterface * solver = model_->solver();
     974    int numberColumns = solver->getNumCols();
     975    // This is number of rows when matrix was passed in
     976    int numberRows = originalNumberRows_;
     977    if (!numberRows)
     978        return 0; // switched off
     979
     980    const double * columnLower = solver->getColLower();
     981    const double * columnUpper = solver->getColUpper();
     982    // modified rhs
     983    double * rhs = CoinCopyOfArray(originalRhs_,numberRows);
     984    // Column copy
     985    const double * element = matrix_.getElements();
     986    const int * row = matrix_.getIndices();
     987    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
     988    const int * columnLength = matrix_.getVectorLengths();
     989    // If bit set then use current
     990    if ((algorithm_&1)!=0) {
     991      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
     992      element = matrix->getElements();
     993      row = matrix->getIndices();
     994      columnStart = matrix->getVectorStarts();
     995      columnLength = matrix->getVectorLengths();
     996      rhs = new double [numberRows];
     997      const double * rowLower = solver->getRowLower();
     998      const double * rowUpper = solver->getRowUpper();
     999      bool good = true;
     1000      for (int iRow = 0; iRow < numberRows; iRow++) {
     1001        if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
     1002          // SOS
     1003          rhs[iRow]=-1.0;
     1004        } else if (rowLower[iRow] > 0.0) {
     1005          good = false;
     1006        } else if (rowUpper[iRow] < 0.0) {
     1007          good = false;
     1008        } else {
     1009          rhs[iRow]=rowUpper[iRow];
     1010        }
     1011      }
     1012      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1013        if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
     1014          good = false;
     1015        CoinBigIndex j;
     1016        int nSOS=0;
     1017        if (!solver->isInteger(iColumn))
     1018          good = false;
     1019        for (j = columnStart[iColumn];
     1020             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1021          if (element[j] < 0.0)
     1022            good = false;
     1023          int iRow = row[j];
     1024          if (rhs[iRow]==-1.0) {
     1025            if (element[j] != 1.0)
     1026              good = false;
     1027            nSOS++;
     1028          }
     1029        }
     1030        if (nSOS!=1)
     1031          good = false;
     1032      }
     1033      if (!good) {
     1034        delete [] rhs;
     1035        setWhen(0); // switch off
     1036        return 0;
     1037      }
     1038    }
     1039    const double * solution = solver->getColSolution();
     1040    const double * objective = solver->getObjCoefficients();
     1041    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1042    double primalTolerance;
     1043    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
     1044
     1045    numRuns_++;
     1046    assert (numberRows == matrix_.getNumRows());
     1047    int iRow, iColumn;
     1048    double direction = solver->getObjSense();
     1049    double * slackCost = new double [numberRows];
     1050    double * modifiedCost = CoinCopyOfArray(objective,numberColumns);
     1051    for (int iRow = 0;iRow < numberRows; iRow++)
     1052      slackCost[iRow]=1.0e30;
     1053    // Take off cost of gub slack
     1054    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1055      if (columnLength[iColumn] == 1) {
     1056        // SOS slack
     1057        double cost = direction*objective[iColumn];
     1058        int iRow = row[columnStart[iColumn]];
     1059        assert (rhs[iRow]<0.0);
     1060        slackCost[iRow]=CoinMin(slackCost[iRow],cost);
     1061      }
     1062    }
     1063    double offset2 = 0.0;
     1064    for (int iRow = 0;iRow < numberRows; iRow++) {
     1065      if( slackCost[iRow] == 1.0e30) {
     1066        slackCost[iRow]=0.0;
     1067      } else {
     1068        offset2 += slackCost[iRow];
     1069        rhs[iRow] = -2.0;
     1070      }
     1071    }
     1072    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1073      double cost = direction * modifiedCost[iColumn];
     1074      CoinBigIndex j;
     1075      for (j = columnStart[iColumn];
     1076           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1077        int iRow = row[j];
     1078        if (rhs[iRow]<0.0) {
     1079          cost -= slackCost[iRow];
     1080        }
     1081      }
     1082      modifiedCost[iColumn] = cost;
     1083    }
     1084    delete [] slackCost;
     1085    double offset;
     1086    solver->getDblParam(OsiObjOffset, offset);
     1087    double newSolutionValue = -offset+offset2;
     1088    int returnCode = 0;
     1089
     1090
     1091    // Get solution array for heuristic solution
     1092    double * newSolution = new double [numberColumns];
     1093    double * rowActivity = new double[numberRows];
     1094    memset(rowActivity, 0, numberRows*sizeof(double));
     1095    if ((algorithm_&2)==0) {
     1096      // get solution as small as possible
     1097      for (iColumn = 0; iColumn < numberColumns; iColumn++)
     1098        newSolution[iColumn] = columnLower[iColumn];
     1099    } else {
     1100      // Get rounded down solution
     1101      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1102        double value = solution[iColumn];
     1103        // Round down integer
     1104        if (fabs(floor(value + 0.5) - value) < integerTolerance) {
     1105          value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn]));
     1106        } else {
     1107          value = CoinMax(floor(value), columnLower[iColumn]);
     1108        }
     1109        // make sure clean
     1110        value = CoinMin(value, columnUpper[iColumn]);
     1111        value = CoinMax(value, columnLower[iColumn]);
     1112        newSolution[iColumn] = value;
     1113      }
     1114    }
     1115    double * contribution = new double [numberColumns];
     1116    int * which = new int [numberColumns];
     1117    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1118      CoinBigIndex j;
     1119      double value = newSolution[iColumn];
     1120      double cost =  modifiedCost[iColumn];
     1121      double forSort = 0.0;
     1122      bool hasSlack=false;
     1123      newSolutionValue += value * cost;
     1124      for (j = columnStart[iColumn];
     1125           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1126        int iRow = row[j];
     1127        rowActivity[iRow] += value * element[j];
     1128        if (rhs[iRow] >0.0) {
     1129          forSort += element[j];
     1130        } else  if (rhs[iRow] == -2.0) {
     1131          hasSlack = true;
     1132        }
     1133      }
     1134      if (forSort && value == 0.0 && columnUpper[iColumn]) {
     1135        if (hasSlack) {
     1136          if (cost>=0.0) {
     1137            forSort = 2.0e30;
     1138          } else {
     1139            forSort = cost/forSort;
     1140          }
     1141        } else {
     1142          forSort = cost/forSort;
     1143        }
     1144      } else {
     1145        // put at end
     1146        forSort = 1.0e30;
     1147      }
     1148      which[iColumn]=iColumn;
     1149      contribution[iColumn]= forSort;
     1150    }
     1151    CoinSort_2(contribution,contribution+numberColumns,which);
     1152    // Go through columns
     1153    for (int jColumn = 0; jColumn < numberColumns; jColumn++) {
     1154      int iColumn = which[jColumn];
     1155      double value = newSolution[iColumn];
     1156      if (value)
     1157        continue;
     1158      bool possible = true;
     1159      CoinBigIndex j;
     1160      for (j = columnStart[iColumn];
     1161           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1162        int iRow = row[j];
     1163        if (rhs[iRow]<0.0&&rowActivity[iRow]) {
     1164          possible = false;
     1165        } else if (rhs[iRow]>=0.0) {
     1166          double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8;
     1167          if (gap<element[j])
     1168            possible = false;
     1169        }
     1170      }
     1171      if (possible) {
     1172        // Increase chosen column
     1173        newSolution[iColumn] = 1.0;
     1174        double cost =  modifiedCost[iColumn];
     1175        newSolutionValue += cost;
     1176        for (CoinBigIndex j = columnStart[iColumn];
     1177             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1178          int iRow = row[j];
     1179          rowActivity[iRow] += element[j];
     1180        }
     1181      }
     1182    }
     1183    if (newSolutionValue < solutionValue) {
     1184        // check feasible
     1185      const double * rowLower = solver->getRowLower();
     1186      const double * rowUpper = solver->getRowUpper();
     1187      memset(rowActivity, 0, numberRows*sizeof(double));
     1188        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1189            CoinBigIndex j;
     1190            double value = newSolution[iColumn];
     1191            if (value) {
     1192                for (j = columnStart[iColumn];
     1193                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1194                    int iRow = row[j];
     1195                    rowActivity[iRow] += value * element[j];
     1196                }
     1197            }
     1198        }
     1199        // check was approximately feasible
     1200        bool feasible = true;
     1201        for (iRow = 0; iRow < numberRows; iRow++) {
     1202            if (rowActivity[iRow] < rowLower[iRow]) {
     1203                if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance)
     1204                    feasible = false;
     1205            } else if (rowActivity[iRow] > rowUpper[iRow]) {
     1206                if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance)
     1207                    feasible = false;
     1208            }
     1209        }
     1210        if (feasible) {
     1211            // new solution
     1212            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
     1213            solutionValue = newSolutionValue;
     1214            //printf("** Solution of %g found by rounding\n",newSolutionValue);
     1215            returnCode = 1;
     1216        } else {
     1217            // Can easily happen
     1218            //printf("Debug CbcHeuristicGreedySOS giving bad solution\n");
     1219        }
     1220    }
     1221    delete [] newSolution;
     1222    delete [] rowActivity;
     1223    delete [] modifiedCost;
     1224    delete [] contribution;
     1225    delete [] which;
     1226    delete [] rhs;
     1227    return returnCode;
     1228}
     1229// update model
     1230void CbcHeuristicGreedySOS::setModel(CbcModel * model)
     1231{
     1232    delete [] originalRhs_;
     1233    gutsOfConstructor(model);
     1234    validate();
     1235}
     1236// Resets stuff if model changes
     1237void
     1238CbcHeuristicGreedySOS::resetModel(CbcModel * model)
     1239{
     1240    delete [] originalRhs_;
     1241    gutsOfConstructor(model);
     1242}
     1243// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
     1244void
     1245CbcHeuristicGreedySOS::validate()
     1246{
     1247    if (model_ && when() < 10) {
     1248        if (model_->numberIntegers() !=
     1249                model_->numberObjects() && (model_->numberObjects() ||
     1250                                            (model_->specialOptions()&1024) == 0)) {
     1251            int numberOdd = 0;
     1252            for (int i = 0; i < model_->numberObjects(); i++) {
     1253                if (!model_->object(i)->canDoHeuristics())
     1254                    numberOdd++;
     1255            }
     1256            if (numberOdd)
     1257                setWhen(0);
     1258        }
     1259        // Only works if coefficients positive and all rows L or SOS
     1260        OsiSolverInterface * solver = model_->solver();
     1261        const double * columnUpper = solver->getColUpper();
     1262        const double * columnLower = solver->getColLower();
     1263        const double * rowLower = solver->getRowLower();
     1264        const double * rowUpper = solver->getRowUpper();
     1265
     1266        int numberRows = solver->getNumRows();
     1267        // Column copy
     1268        const double * element = matrix_.getElements();
     1269        const int * row = matrix_.getIndices();
     1270        const CoinBigIndex * columnStart = matrix_.getVectorStarts();
     1271        const int * columnLength = matrix_.getVectorLengths();
     1272        bool good = true;
     1273        assert (originalRhs_);
     1274        for (int iRow = 0; iRow < numberRows; iRow++) {
     1275          if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) {
     1276            // SOS
     1277            originalRhs_[iRow]=-1.0;
     1278          } else if (rowLower[iRow] > 0.0) {
     1279                good = false;
     1280          } else if (rowUpper[iRow] < 0.0) {
     1281                good = false;
     1282          } else {
     1283            originalRhs_[iRow]=rowUpper[iRow];
     1284          }
     1285        }
     1286        int numberColumns = solver->getNumCols();
     1287        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1288            if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0)
     1289                good = false;
     1290            CoinBigIndex j;
     1291            int nSOS=0;
     1292            if (!solver->isInteger(iColumn))
     1293              good = false;
     1294            for (j = columnStart[iColumn];
     1295                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     1296                if (element[j] < 0.0)
     1297                    good = false;
     1298                int iRow = row[j];
     1299                if (originalRhs_[iRow]==-1.0) {
     1300                  if (element[j] != 1.0)
     1301                    good = false;
     1302                  nSOS++;
     1303                }
     1304            }
     1305            if (nSOS!=1)
     1306              good = false;
     1307        }
     1308        if (!good)
     1309            setWhen(0); // switch off
     1310    }
     1311}
Note: See TracChangeset for help on using the changeset viewer.