Order path separator

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<std::pair<HighsInt, double>> 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<std::pair<double, double>> 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<HighsInt, double>& i,
+                const std::pair<HighsInt, double>& j) {
+              return rowscore[i.first] > rowscore[j.first];
+            });
+    pdqsort(outArcRows.begin() + colOutArcs[col].first,
+            outArcRows.begin() + colOutArcs[col].second,
+            [&](const std::pair<HighsInt, double>& i,
+                const std::pair<HighsInt, double>& j) {
+              return rowscore[i.first] > rowscore[j.first];
+            });
   }
 
   HighsCutGeneration cutGen(lpRelaxation, cutpool);
@@ -258,31 +295,10 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation,
               const std::vector<std::pair<HighsInt, HighsInt>>& colArcs,
               const std::vector<std::pair<HighsInt, double>>& 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;

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];
 

highs/mip/HighsTransformedLp.h

@@ -37,6 +37,7 @@ class HighsTransformedLp {
   std::vector<double> lbDist;
   std::vector<double> ubDist;
   std::vector<double> boundDist;
+  std::vector<double> 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<double>& vals, std::vector<double>& upper,
                  std::vector<double>& solval, std::vector<HighsInt>& inds,
                  double& rhs, bool& integralPositive, bool preferVbds = false);