Changeset 1564 for trunk


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

add some heuristic variants

Location:
trunk/Cbc/src
Files:
4 edited

Legend:

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

    r1505 r1564  
    249249    inline void setWhereFrom(int value) {
    250250        whereFrom_ = value;
     251    }
     252    inline int whereFrom() const {
     253        return whereFrom_;
    251254    }
    252255    /** Upto this depth we call the tree shallow and the heuristic can be called
  • 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}
  • trunk/Cbc/src/CbcHeuristicGreedy.hpp

    r1432 r1564  
    186186};
    187187
     188/** Greedy heuristic for SOS and L rows (and positive elements)
     189 */
     190
     191class CbcHeuristicGreedySOS : public CbcHeuristic {
     192public:
     193
     194    // Default Constructor
     195    CbcHeuristicGreedySOS ();
     196
     197    /* Constructor with model - assumed before cuts
     198       Initial version does not do Lps
     199    */
     200    CbcHeuristicGreedySOS (CbcModel & model);
     201
     202    // Copy constructor
     203    CbcHeuristicGreedySOS ( const CbcHeuristicGreedySOS &);
     204
     205    // Destructor
     206    ~CbcHeuristicGreedySOS ();
     207
     208    /// Clone
     209    virtual CbcHeuristic * clone() const;
     210    /// Assignment operator
     211    CbcHeuristicGreedySOS & operator=(const CbcHeuristicGreedySOS& rhs);
     212    /// Create C++ lines to get to current state
     213    virtual void generateCpp( FILE * fp) ;
     214
     215    /// update model (This is needed if cliques update matrix etc)
     216    virtual void setModel(CbcModel * model);
     217
     218    using CbcHeuristic::solution ;
     219    /** returns 0 if no solution, 1 if valid solution.
     220        Sets solution values if good, sets objective value (only if good)
     221        We leave all variables which are at one at this node of the
     222        tree to that value and will
     223        initially set all others to zero.  We then sort all variables in order of their cost
     224        divided by the number of entries in rows which are not yet covered.  We randomize that
     225        value a bit so that ties will be broken in different ways on different runs of the heuristic.
     226        We then choose the best one and set it to one and repeat the exercise.
     227
     228    */
     229    virtual int solution(double & objectiveValue,
     230                         double * newSolution);
     231    /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
     232    virtual void validate() ;
     233    /// Resets stuff if model changes
     234    virtual void resetModel(CbcModel * model);
     235    /* Algorithm
     236       Bits
     237       1 bit - use current model, otherwise original
     238       2 - use current solution as starting point, otherwise pure greedy
     239       4 - use duals to modify greedy
     240       8 - use duals on GUB/SOS in special way
     241    */
     242    inline int algorithm() const {
     243        return algorithm_;
     244    }
     245    inline void setAlgorithm(int value) {
     246        algorithm_ = value;
     247    }
     248    // Only do this many times
     249    inline int numberTimes() const {
     250        return numberTimes_;
     251    }
     252    inline void setNumberTimes(int value) {
     253        numberTimes_ = value;
     254    }
     255
     256protected:
     257    /// Guts of constructor from a CbcModel
     258    void gutsOfConstructor(CbcModel * model);
     259    // Data
     260
     261    // Original RHS - if -1.0 then SOS otherwise <= value
     262    double * originalRhs_;
     263    // Original matrix by column
     264    CoinPackedMatrix matrix_;
     265    // original number of rows
     266    int originalNumberRows_;
     267    /* Algorithm
     268    */
     269    int algorithm_;
     270    /// Do this many times
     271    int numberTimes_;
     272
     273};
     274
    188275
    189276#endif
  • trunk/Cbc/src/CbcHeuristicRENS.cpp

    r1499 r1564  
    1414#include "CbcHeuristicRENS.hpp"
    1515#include "CoinWarmStartBasis.hpp"
     16#include "CoinSort.hpp"
    1617#include "CbcBranchActual.hpp"
    1718#include "CbcStrategy.hpp"
     
    8990
    9091    const double * currentSolution = solver->getColSolution();
    91     const double * dj = solver->getReducedCost();
     92    int numberColumns = solver->getNumCols();
     93    double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
    9294    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
    93     int numberColumns = newSolver->getNumCols();
    9495    double direction = newSolver->getObjSense();
    9596    int type = rensType_&15;
     
    9798    const double * colLower = newSolver->getColLower();
    9899    const double * colUpper = newSolver->getColUpper();
    99     if ((type&3)==3) {
    100         double total=0.0;
    101         int n=0;
    102         CoinWarmStartBasis * basis =
    103           dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
    104         if (basis) {
    105             for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    106               if (colUpper[iColumn]>colLower[iColumn]&&
    107                   basis->getStructStatus(iColumn) !=
    108                   CoinWarmStartBasis::basic) {
    109                   n++;
    110                   total += fabs(dj[iColumn]);
     100    int numberFixed = 0;
     101    if ((type&7)==3) {
     102      double total=0.0;
     103      int n=0;
     104      CoinWarmStartBasis * basis =
     105        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     106      if (basis&&basis->getNumArtificial()) {
     107        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     108          if (colUpper[iColumn]>colLower[iColumn]&&
     109              basis->getStructStatus(iColumn) !=
     110              CoinWarmStartBasis::basic) {
     111            n++;
     112            total += fabs(dj[iColumn]);
     113          }
     114        }
     115        if (n)
     116          djTolerance = (0.01*total)/static_cast<double>(n);
     117        delete basis;
     118      }
     119    } else if ((type&7)==4) {
     120      double * sort = new double [numberColumns];
     121      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     122        sort[iColumn]=1.0e30;
     123        if (colUpper[iColumn]>colLower[iColumn]) {
     124          sort[iColumn] = fabs(dj[iColumn]);
     125        }
     126      }
     127      std::sort(sort,sort+numberColumns);
     128      int last = static_cast<int>(numberColumns*fractionSmall_);
     129      djTolerance = sort[last];
     130      delete [] sort;
     131    } else if ((type&7)==5||(type&7)==6) {
     132      // SOS type fixing
     133      bool fixSets = (type&7)==5;
     134      CoinWarmStartBasis * basis =
     135        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
     136      if (basis&&basis->getNumArtificial()) {
     137        //const double * rowLower = solver->getRowLower();
     138        const double * rowUpper = solver->getRowUpper();
     139       
     140        int numberRows = solver->getNumRows();
     141        // Column copy
     142        const CoinPackedMatrix * matrix = solver->getMatrixByCol();
     143        const double * element = matrix->getElements();
     144        const int * row = matrix->getIndices();
     145        const CoinBigIndex * columnStart = matrix->getVectorStarts();
     146        const int * columnLength = matrix->getVectorLengths();
     147        double * bestDj = new double [numberRows];
     148        for (int i=0;i<numberRows;i++) {
     149          if (rowUpper[i]==1.0)
     150            bestDj[i]=1.0e20;
     151          else
     152            bestDj[i]=1.0e30;
     153        }
     154        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     155          if (colUpper[iColumn]>colLower[iColumn]) {
     156            CoinBigIndex j;
     157            if (currentSolution[iColumn]>1.0e-6&&
     158                currentSolution[iColumn]<0.999999) {
     159              for (j = columnStart[iColumn];
     160                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     161                int iRow = row[j];
     162                bestDj[iRow]=1.0e30;
     163              }
     164            } else if ( basis->getStructStatus(iColumn) !=
     165              CoinWarmStartBasis::basic) {
     166              for (j = columnStart[iColumn];
     167                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     168                int iRow = row[j];
     169                if (bestDj[iRow]<1.0e30) {
     170                  if (element[j] != 1.0)
     171                    bestDj[iRow]=1.0e30;
     172                  else
     173                    bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]);
     174                }
    111175              }
    112176            }
    113             if (n)
    114                 djTolerance = (0.01*total)/static_cast<double>(n);
    115             delete basis;
    116         }
     177          }
     178        }
     179        int nSOS=0;
     180        double * sort = new double [numberRows];
     181        for (int i=0;i<numberRows;i++) {
     182          if (bestDj[i]<1.0e30) {
     183            sort[nSOS++]=bestDj[i];
     184          }
     185        }
     186        if (10*nSOS>8*numberRows) {
     187          std::sort(sort,sort+nSOS);
     188          int last = static_cast<int>(nSOS*0.9*fractionSmall_);
     189          double tolerance = sort[last];
     190          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     191            if (colUpper[iColumn]>colLower[iColumn]) {
     192              CoinBigIndex j;
     193              if (currentSolution[iColumn]<=1.0e-6||
     194                  currentSolution[iColumn]>=0.999999) {
     195                if (fixSets) {
     196                  for (j = columnStart[iColumn];
     197                       j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     198                    int iRow = row[j];
     199                    if (bestDj[iRow]<1.0e30&&bestDj[iRow]>=tolerance) {
     200                      numberFixed++;
     201                      if (currentSolution[iColumn]<=1.0e-6)
     202                        newSolver->setColUpper(iColumn,0.0);
     203                      else if (currentSolution[iColumn]>=0.999999)
     204                        newSolver->setColLower(iColumn,1.0);
     205                    }
     206                  }
     207                } else if (columnLength[iColumn]==1) {
     208                  // leave more slacks
     209                  int iRow = row[columnStart[iColumn]];
     210                  if (bestDj[iRow]<1.0e30) {
     211                    // fake dj
     212                    dj[iColumn] *= 0.000001;
     213                  }
     214                }
     215              }
     216            }
     217          }
     218          if (fixSets)
     219            djTolerance = 1.0e30;
     220          else
     221            djTolerance = tolerance;
     222        }
     223        delete basis;
     224        delete [] sort;
     225        delete [] bestDj;
     226      }
    117227    }
    118 
     228   
    119229    double primalTolerance;
    120230    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
    121231
    122232    int i;
    123     int numberFixed = 0;
    124233    int numberTightened = 0;
    125234    int numberAtBound = 0;
     
    183292#endif
    184293    }
     294    delete [] dj;
    185295    if (numberFixed > numberIntegers / 5) {
    186296        if ( numberFixed < numberColumns / 5) {
Note: See TracChangeset for help on using the changeset viewer.