Finish draft

Opt-Mucca · Opt-Mucca · commit 9d4b2e5a406b · 2025-11-18T12:21:08.000+01:00
diff --git a/highs/mip/HighsPathSeparator.cpp b/highs/mip/HighsPathSeparator.cpp
@@ -15,6 +15,7 @@
 #include "mip/HighsLpRelaxation.h"
 #include "mip/HighsMipSolverData.h"
 #include "mip/HighsTransformedLp.h"
+#include "pdqsort.h"
 
 enum class RowType : int8_t {
   kUnusuable = -2,
@@ -58,10 +59,7 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation,
   }
 
   std::vector<HighsInt> numContinuous(lp.num_row_);
-  // Score will be used for deciding order in which rows get aggregated
-  // TODO: The second entry is currently a hacky way to store the norm
-  std::vector<std::pair<double, double>> rowscore(lp.num_row_,
-                                                  std::make_pair(0.0, 0.0));
+
   size_t maxAggrRowSize = 0;
   for (HighsInt col : mip.mipdata_->continuous_cols) {
     if (transLp.boundDistance(col) == 0.0) continue;
@@ -70,10 +68,6 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation,
     for (HighsInt i = lp.a_matrix_.start_[col];
          i != lp.a_matrix_.start_[col + 1]; ++i) {
       ++numContinuous[lp.a_matrix_.index_[i]];
-      // Add the fractional score of the row
-      rowscore[i].first +=
-          lp.a_matrix_.value_[i] * transLp.getFracVbEstimate(col);
-      rowscore[i].second += lp.a_matrix_.value_[i];
     }
   }
 
@@ -114,6 +108,25 @@ 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];
+      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 > 0) {
+      rowscore[row].first /= rowscore[row].second;
+      rowscore[row].second = 1;
+    }
+  }
+
   // 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
@@ -170,6 +183,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);
@@ -266,31 +293,9 @@ 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;
diff --git 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] <= simpleLbDist[col] &&
+        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] <= simpleUbDist[col] &&
+               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
@@ -16,7 +16,6 @@
 
 #include <vector>
 
-#include "HighsLpRelaxation.h"
 #include "lp_data/HConst.h"
 #include "mip/HighsImplications.h"
 #include "util/HighsCDouble.h"
@@ -38,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,
@@ -53,28 +53,8 @@ class HighsTransformedLp {
 
   double boundDistance(HighsInt col) const { return boundDist[col]; }
 
-  // TODO: Should this simple be calculated when we initialise the transLp?
-  double getFracVbEstimate(HighsInt col) const {
-    HighsInt vbCol = -1;
-    double vbCoef = 0.0;
-    if (ubDist[col] <= lbDist[col] && bestVlb[col].first != -1) {
-      vbCol = bestVlb[col].first;
-      vbCoef = bestVlb[col].second.coef;
-    } else if (lbDist[col] <= ubDist[col] && bestVub[col].first != -1) {
-      vbCol = bestVub[col].first;
-      vbCoef = bestVub[col].second.coef;
-    } else if (bestVlb[col].first != -1) {
-      vbCol = bestVlb[col].first;
-      vbCoef = bestVlb[col].second.coef;
-    } else if (bestVub[col].first != -1) {
-      vbCol = bestVub[col].first;
-      vbCoef = bestVub[col].second.coef;
-    }
-    if (vbCol == -1) return 0;
-    // TODO: This differs from the SCIP definition? Their frac can be > 1?
-    double val = vbCoef * lprelaxation.solutionValue(col);
-    double frac = std::max(val - std::floor(val), 0.0);
-    return std::min(frac, 1 - frac);
+  double getColFractionality(HighsInt col) const {
+    return colFractionality[col];
   }
 
   bool transform(std::vector<double>& vals, std::vector<double>& upper,
