diff --git a/highs/mip/HighsPathSeparator.cpp b/highs/mip/HighsPathSeparator.cpp index f3b6df3062..f3bdf4d5c4 100644 --- a/highs/mip/HighsPathSeparator.cpp +++ b/highs/mip/HighsPathSeparator.cpp @@ -10,6 +10,7 @@ #include "mip/HighsPathSeparator.h" +#include "../extern/pdqsort/pdqsort.h" #include "mip/HighsCutGeneration.h" #include "mip/HighsLpAggregator.h" #include "mip/HighsLpRelaxation.h" @@ -65,8 +66,9 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, maxAggrRowSize += lp.a_matrix_.start_[col + 1] - lp.a_matrix_.start_[col]; for (HighsInt i = lp.a_matrix_.start_[col]; - i != lp.a_matrix_.start_[col + 1]; ++i) + i != lp.a_matrix_.start_[col + 1]; ++i) { ++numContinuous[lp.a_matrix_.index_[i]]; + } } std::vector> colSubstitutions( @@ -106,6 +108,27 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, rowtype[i] = RowType::kUnusuable; } + // Score will be used for deciding order in which rows get aggregated + std::vector> rowscore(lp.num_row_, + std::make_pair(0.0, 0.0)); + for (HighsInt col = 0; col != lp.num_col_; ++col) { + for (HighsInt i = lp.a_matrix_.start_[col]; + i != lp.a_matrix_.start_[col + 1]; ++i) { + HighsInt row = lp.a_matrix_.index_[i]; + if (rowtype[row] == RowType::kUnusuable) continue; + double val = lp.a_matrix_.value_[i]; + rowscore[row].first += std::abs(val * transLp.getColFractionality(col)); + rowscore[row].second += std::abs(val); + } + } + for (HighsInt row = 0; row != lp.num_row_; ++row) { + if (rowscore[row].second > mip.mipdata_->feastol) { + rowscore[row].first /= rowscore[row].second; + } else { + rowscore[row].first = 0.0; + } + } + // for each continuous variable with nonzero transformed solution value // remember the <= and == rows where it is present with a positive coefficient // in its set of in-arc rows. Treat >= rows as <= rows with reversed sign @@ -162,6 +185,20 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, colInArcs[col].second = inArcRows.size(); colOutArcs[col].second = outArcRows.size(); + + // Sort the in and outgoing arcs by their scores + pdqsort(inArcRows.begin() + colInArcs[col].first, + inArcRows.begin() + colInArcs[col].second, + [&](const std::pair& i, + const std::pair& j) { + return rowscore[i.first] > rowscore[j.first]; + }); + pdqsort(outArcRows.begin() + colOutArcs[col].first, + outArcRows.begin() + colOutArcs[col].second, + [&](const std::pair& i, + const std::pair& j) { + return rowscore[i.first] > rowscore[j.first]; + }); } HighsCutGeneration cutGen(lpRelaxation, cutpool); @@ -258,31 +295,10 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, const std::vector>& colArcs, const std::vector>& arcRows, HighsInt& row, double& weight) { - HighsInt arcRow = randgen.integer(colArcs[bestArcCol].first, - colArcs[bestArcCol].second); - HighsInt r = arcRows[arcRow].first; - double w = -val / arcRows[arcRow].second; - if (!isRowInCurrentPath(r) && checkWeight(w)) { - row = r; - weight = w; - return true; - } - - for (HighsInt nextRow = arcRow + 1; - nextRow < colArcs[bestArcCol].second; ++nextRow) { - r = arcRows[nextRow].first; - w = -val / arcRows[nextRow].second; - if (!isRowInCurrentPath(r) && checkWeight(w)) { - row = r; - weight = w; - return true; - } - } - - for (HighsInt nextRow = colArcs[bestArcCol].first; nextRow < arcRow; - ++nextRow) { - r = arcRows[nextRow].first; - w = -val / arcRows[nextRow].second; + for (HighsInt arcRow = colArcs[bestArcCol].first; + arcRow != colArcs[bestArcCol].second; ++arcRow) { + HighsInt r = arcRows[arcRow].first; + double w = -val / arcRows[arcRow].second; if (!isRowInCurrentPath(r) && checkWeight(w)) { row = r; weight = w; diff --git a/highs/mip/HighsTransformedLp.cpp b/highs/mip/HighsTransformedLp.cpp index 703b377e09..6050f0b816 100644 --- a/highs/mip/HighsTransformedLp.cpp +++ b/highs/mip/HighsTransformedLp.cpp @@ -32,6 +32,7 @@ HighsTransformedLp::HighsTransformedLp(const HighsLpRelaxation& lprelaxation, std::make_pair(-1, HighsImplications::VarBound())); boundTypes.resize(numTransformedCol); vectorsum.setDimension(numTransformedCol); + colFractionality.resize(lprelaxation.numCols()); for (HighsInt col : mipsolver.mipdata_->continuous_cols) { mipsolver.mipdata_->implications.cleanupVarbounds(col); @@ -57,9 +58,26 @@ HighsTransformedLp::HighsTransformedLp(const HighsLpRelaxation& lprelaxation, if (ubDist[col] <= mipsolver.mipdata_->feastol) ubDist[col] = 0.0; boundDist[col] = std::min(lbDist[col], ubDist[col]); + + if (lbDist[col] <= ubDist[col] && lbDist[col] == 0 && + bestVlb[col].first != -1) { + const double frac = + std::min(lpSolution.col_value[bestVlb[col].first], + 1 - lpSolution.col_value[bestVlb[col].first]); + colFractionality[col] = bestVlb[col].second.coef * frac; + } else if (ubDist[col] <= lbDist[col] && ubDist[col] == 0 && + bestVub[col].first != -1) { + const double frac = + std::min(lpSolution.col_value[bestVub[col].first], + 1 - lpSolution.col_value[bestVub[col].first]); + colFractionality[col] = bestVub[col].second.coef * frac; + } } for (HighsInt col : mipsolver.mipdata_->integral_cols) { + double frac = + lpSolution.col_value[col] - std::floor(lpSolution.col_value[col]); + colFractionality[col] = std::min(frac, 1 - frac); double bestub = mipsolver.mipdata_->domain.col_upper_[col]; double bestlb = mipsolver.mipdata_->domain.col_lower_[col]; diff --git a/highs/mip/HighsTransformedLp.h b/highs/mip/HighsTransformedLp.h index 5c3d2d01db..7a49aaff84 100644 --- a/highs/mip/HighsTransformedLp.h +++ b/highs/mip/HighsTransformedLp.h @@ -37,6 +37,7 @@ class HighsTransformedLp { std::vector lbDist; std::vector ubDist; std::vector boundDist; + std::vector colFractionality; enum class BoundType : uint8_t { kSimpleUb, kSimpleLb, @@ -52,6 +53,10 @@ class HighsTransformedLp { double boundDistance(HighsInt col) const { return boundDist[col]; } + double getColFractionality(HighsInt col) const { + return colFractionality[col]; + } + bool transform(std::vector& vals, std::vector& upper, std::vector& solval, std::vector& inds, double& rhs, bool& integralPositive, bool preferVbds = false);